As France enters its second heatwave of 2026, can we produce more detailed plots than the excellent visualizations provided by ShowYourStripes?

MétéoFrance offers its monthly SIM2 dataset albeit over a shorter time span (currently 1970–2025). The dataset includes temperature, precipitation and other variables on an 8 km resolution grid.

We will select a French city, retrieve its geographic coordinates, build the grid for a specific month over the 1970–2025 period, extract the data from the grid at that location and plot the temperature anomaly.

We will use {terra} to create the grid from the tabular files containing cell centers and weather variables, and terra::extract() to get all temperatures.

library(tidyverse)
library(fs)
library(janitor)
library(osmdata)
library(sf)
library(terra)
library(glue)
library(memoise)

invisible(
  Sys.setlocale(category = "LC_ALL", 
                locale = "en_GB.UTF8"))


#' Geolocate using OSM (Nominatim API)
#'
#' using first result
#' Memoized function
#'
#' @param location (char): place name (geocodable via OSM)
#'
#' @returns (SpatVector): first result
geolocate <- memoise(function(location) {
  loc <- getbb(location, format_out = "sf_polygon")  |> 
    st_point_on_surface() |>
    slice_head(n = 1)
  
  message(glue("{loc$display_name} — WGS84 : {loc$geometry}"))
  
  return(loc|>
           st_transform("IGNF:NTFLAMB2E") |>
           as("SpatVector"))
})

#' Generate a monthly temperature chart since 1970
#'
#' @param sim2 (data.frame): Météo-France SIM2 data over the period
#' @param month (char): month number "01"..."12"
#' @param location (char): place name (geocodable via OSM); memoized
#' @param output_dir (char): directory path where a PNG file will be written, if not NULL
#'
#' @returns (ggplot and optionally a file on disk)
generate_chart <- function(
    sim2,
    month,
    location,
    output_dir = NULL) {
  stopifnot(month %in% sprintf("%02d", 1:12))
  month_name <- format(ymd(glue("0000-{month}-01")), "%B")
  
  sim2_raster <- sim2 |>
    filter(str_detect(date, glue("{month}$"))) |>
    mutate(
      x = lambx * 100,
      y = lamby * 100,
      layer = date,
      temp = t,
      .keep = "none") |>
    rast(
      type = "xylz",
      crs = "IGNF:NTFLAMB2E")
  
  loc <- geolocate(location)
  
  temperatures <- sim2_raster |>
    terra::extract(loc) |>
    select(-ID) |>
    pivot_longer(
      cols = everything(),
      names_to = "month",
      values_to = "temperature") |>
    mutate(
      year = as.integer(str_sub(month, 1, 4)),
      anomaly = temperature - mean(temperature[year >= 1991 & year <= 2020],
                                   na.rm = TRUE))
  
  p <- temperatures |>
    ggplot(aes(year, anomaly)) +
    geom_col(aes(fill = anomaly)) +
    geom_smooth(method = "loess",
                formula = y ~ x) +
    scale_fill_gradient2(
      high = scales::muted("red"),
      mid = "white",
      low = scales::muted("blue")) +
    scale_x_continuous(breaks = scales::breaks_pretty()) +
    scale_y_continuous(breaks = scales::breaks_pretty()) +
    labs(
      title = glue("Average monthly anomaly temperature — {month_name}"),
      subtitle = location,
      x = "year",
      y = "departure from average* (°C)",
      fill = "°C",
      caption = glue(
        "https://r.iresmi.net/ — {Sys.Date()}
         data: Météo-France SIM2 — *baseline: 1991–2020 normal for {month_name}")) +
    theme(
      text = element_text(family = "Ubuntu"),
      plot.caption = element_text(size = 7))
  
  if (!is.null(output_dir)) {
    dir_create(output_dir)
    
    ggsave(
      glue("{output_dir}/tm_{month}_{make_clean_names(location)}.png"),
      plot = p,
      width = 20,
      height = 20 / 1.618,
      units = "cm",
      dpi = 150)
  }
  
  return(p)
}

The data is a bunch of compressed CSV.

# https://meteo.data.gouv.fr/datasets/65e040c50a5c6872ebebc711
# Climate change data - monthly SIM
# all files MENS_SIM2_*-*.csv.gz
sim2 <- dir_ls("data") |>
  read_delim(
    delim = ";",
    locale = locale(decimal_mark = "."),
    name_repair = make_clean_names)

Now we just call our function.

generate_chart(sim2,
               month = "06",
               location = "Paris, France")

If we want each month:

# for each month
sprintf("%02d", 1:12) |>
  map(\(x) generate_chart(sim2,
                           month = x,
                           location = "Grenoble, France",
                           output_dir = "results"), 
  .progress = TRUE)

Note that scales are not constant across plots; if we want to compare month (or places) we should fix the y-axis and the color scale. It’s left as an exercise to the reader if you want to make a nice poster…

Join our workshop on Introduction to Bayesian Multiple Imputation with the rblimp package,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: Introduction to Bayesian Multiple Imputation with the rblimp package

Date: Thursday, July 23rd, 18:00 – 20:00 CEST (Rome, Berlin, Paris timezone) 

Speaker: Ermioni Athanasiadi is a PhD student at the University of Siegen, holding a master’s degree in psychology from University of Tübingen and currently pursuing a Master’s degree in Statistics and Data Science at Hasselt University, Belgium. Her research focuses on missing data methods in small-sample settings.  She is also passionate about interdisciplinary perspectives on statistics and research methodology.

Description:   Multiple Imputation has become the gold standard for handling missing data in applied research, with many researchers making use of the mice package. Blimp is an alternative, flexible software that allows multiple imputation within a fully Bayesian framework. In this workshop, we will cover the basics of multiple imputation and Blimp’s modeling framework. We will then explore how to specify imputation models, incorporate auxiliary variables, assess convergence and post-imputation diagnostics and conduct pooled analyses using multiply imputed datasets.

The workshop may be helpful both for those who are new to multiple imputation, and for those who have previously used multiple imputation methods and would like to learn about a Bayesian alternative.

Participants should install the Blimp software beforehand: https://www.appliedmissingdata.com/blimp

Minimal registration fee: 20 euro (or 20 USD or 800 UAH)

Please note that the registration confirmation is sent 1 day before the workshop to all registered participants rather than immediately after registration

How can I register?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or  https://bit.ly/3PFxtNA and donate at least 20 euro. Feel free to donate more if you can, all proceeds go directly to support Ukraine.

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the registration form, attaching a screenshot of a donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after donation).

If you are not personally interested in attending, you can also contribute by sponsoring a participation of a student, who will then be able to participate for free. If you choose to sponsor a student, all proceeds will also go directly to organisations working in Ukraine. You can either sponsor a particular student or you can leave it up to us so that we can allocate the sponsored place to students who have signed up for the waiting list.

How can I sponsor a student?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or https://bit.ly/3PFxtNA and donate at least 20 euro (or 17 GBP or 20 USD or 800 UAH). Feel free to donate more if you can, all proceeds go to support Ukraine!

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the sponsorship form, attaching the screenshot of the donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after the donation). You can indicate whether you want to sponsor a particular student or we can allocate this spot ourselves to the students from the waiting list. You can also indicate whether you prefer us to prioritize students from developing countries when assigning place(s) that you sponsored.

If you are a university student and cannot afford the registration fee, you can also sign up for the waiting list here. (Note that you are not guaranteed to participate by signing up for the waiting list).

You can also find more information about this workshop series,  a schedule of our future workshops as well as a list of our past workshops which you can get the recordings & materials here.

Looking forward to seeing you during the workshop!

Flowcharts should be beautiful. Just like this CC photo from Wasif Malik,

Thanks to Alan Haynes and
his excellent suggestions, I have spent some time improving the
flowchart component of the Gmisc package. The result is not meant to be
another decorative diagram tool. It is meant for the kind of figures
researchers keep redrawing by hand: CONSORT diagrams, cohort derivation
charts, screening flows, data-cleaning audit trails, and the small but
important maps that explain how a study population came to be.

I like tools such as Excalidraw
for thinking. They are fast, expressive, and excellent for
conversations. But when a figure enters a manuscript, the needs change.
Counts must be updated. Exclusions must match the analysis script.
Treatment arms should align. Follow-up losses should be traceable. The
figure should survive reviewer round three without becoming a manual
editing project.

That is the space where flowchart() in Gmisc is useful:
the diagram becomes part of the research workflow.

The figure above is the kind of chart I want Gmisc to make feel
natural. It is still a grid graphic in R, but it has the visual grammar
of a manuscript figure: grouped arms, side exclusions, count badges,
phase labels, and arrows that do not need nudging after every text
change.

Every figure in this post is generated by code, and the code is
included below each image. They all share the same two-line
preamble:

library(Gmisc)
library(grid)

To save any of them to a file, wrap the call in a graphics device,
e.g.

png("01-consort-color.png", width = 9, height = 7, units = "in", res = 180, bg = "white")
# ... the flowchart code ...
dev.off()

The CONSORT figure above is produced by:

options(boxGrobTxtPadding = unit(3, "mm"))

box_fill <- gpar(fill = "#DDEEFF", col = "#336699", lwd = 1.5)
con_gp <- gpar(col = "#336699", lwd = 1.5, fill = "#336699")
side_gp <- gpar(col = "#CC8800", lwd = 1.2, fill = "#CC8800")
excl_fill <- gpar(fill = "#FFF8E1", col = "#CC8800", lwd = 1.2)
heading_gp <- gpar(fill = "#C8DAF7", col = "#2F5F9F", lwd = 1.1)
badge_gp <- gpar(fill = "#336699", col = NA)
badge_txt_gp <- gpar(col = "white", cex = 0.65)

main_arm_margin <- 0.28
main_x <- 0.5
exclusion_margin <- 0.05

grid.newpage()
flowchart(
  assessed = boxGrob(
    "Patients assessed for eligibility",
    x = main_x, box_gp = box_fill,
    badge_label = "840", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
  ),
  randomised = boxGrob(
    "Randomised",
    x = main_x, box_gp = box_fill,
    badge_label = "126", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
  ),
  arms = list(
    cast = boxGrob(
      "Randomised to\ncast immobilisation",
      box_gp = box_fill,
      badge_label = "62", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
    ),
    surgical = boxGrob(
      "Randomised to\nsurgery",
      box_gp = box_fill,
      badge_label = "64", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
    )
  ),
  lost = list(
    lost_cast = boxGrob(
      "Lost to follow-up (n = 2)\n  1 no response\n  1 other surgery",
      just = "left", box_gp = excl_fill
    ),
    lost_surgical = boxGrob(
      "Lost to follow-up (n = 3)\n  2 no response\n  1 other surgery",
      just = "left", box_gp = excl_fill
    )
  ),
  analysis = list(
    analysis_cast = boxGrob(
      "Included in\nprimary analysis",
      box_gp = box_fill,
      badge_label = "60", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
    ),
    analysis_surgical = boxGrob(
      "Included in\nprimary analysis",
      box_gp = box_fill,
      badge_label = "61", badge_gp = badge_gp, badge_txt_gp = badge_txt_gp
    )
  )
) |>
  spread(axis = "y", margin = unit(5, "mm"), exclude = "lost") |>
  align(
    axis = "y",
    subelement = "lost",
    references = list("arms", "analysis")
  ) |>
  equalizeWidths(subelement = list("arms", "analysis")) |>
  spread(axis = "x", subelement = "arms", margin = main_arm_margin) |>
  spread(axis = "x", subelement = "analysis", margin = main_arm_margin) |>
  spread(axis = "x", subelement = "lost", margin = exclusion_margin) |>
  phaseLabel("arms", "Allocation", box_gp = heading_gp) |>
  phaseLabel("analysis", "Analysis", box_gp = heading_gp) |>
  insert(list(excluded = boxGrob(
    "Excluded (n = 714)\n  477 stable ankle mortise\n   64 incongruent ankle mortise\n   30 previous serious trauma\n  143 other reasons",
    just = "left", box_gp = excl_fill
  )), after = "assessed") |>
  move(subelement = "excluded", x = 1 - exclusion_margin, just = "right") |>
  align(
    axis = "y",
    subelement = "excluded",
    references = list("assessed", "randomised")
  ) |>
  connect("assessed", "excluded", type = "L", lty_gp = side_gp, arrow_size = 3, smooth = TRUE) |>
  connect("randomised", "arms", type = "N", lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  connect("assessed", "randomised", type = "v", lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  connect("arms", "lost", type = "L", lty_gp = side_gp, arrow_size = 3, smooth = TRUE) |>
  connect("arms", "analysis", type = "v", lty_gp = con_gp, arrow_size = 3) |>
  print()

A figure that can
change with the analysis

The biggest advantage of drawing a flowchart in code is not that code
is elegant. It is that research figures are rarely finished when we
think they are.

The inclusion count changes after a database refresh. A reviewer asks
for a sensitivity analysis. Someone notices that two exclusion
categories should be split. The statistician reruns the cohort
definition. If the diagram is hand-drawn, every one of those changes
creates a small risk of mismatch between the paper and the actual
analysis.

If the chart is generated, it can sit beside the code that produced
the numbers.

flowchart(...) |>
  spread(axis = "y") |>
  spread(subelement = "arms", axis = "x") |>
  connect("randomised", "arms", type = "N")

That is the mental model: define boxes, arrange boxes, connect boxes.
The final result can still be polished, but it remains reproducible.

Cohort
derivation from data people already have

Most clinical researchers do not start with a perfect trial flow.
They start with a registry extract, an EHR table, a REDCap project, an
Excel sheet from a collaborator, or a combination of all of them.

That workflow deserves a clear figure too.

This kind of diagram is useful because it does not only show who was
included. It shows how the study base was assembled: what sources were
linked, where exclusions entered, and which analytic populations came
out at the end.

I find this especially helpful for observational studies. A table can
report baseline characteristics, but a flowchart explains the
construction of the cohort. It gives the reader a quick answer to: “What
happened between the raw data and the model?”

source_gp <- gpar(fill = "#E8F5E9", col = "#2E7D32", lwd = 1.4)
link_gp <- gpar(fill = "#E3F2FD", col = "#1565C0", lwd = 1.4)
cohort_gp <- gpar(fill = "#FFF8E1", col = "#C69214", lwd = 1.4)
side_gp <- gpar(fill = "#FCE4EC", col = "#AD1457", lwd = 1.2)
final_gp <- gpar(fill = "#EDE7F6", col = "#512DA8", lwd = 1.4)
con_gp <- gpar(col = "#455A64", fill = "#455A64", lwd = 1.4)
excl_gp <- gpar(col = "#AD1457", fill = "#AD1457", lwd = 1.2)

source_margin <- 0.05
output_margin <- 0.05
main_x <- 0.5
main_path <- list("linked", "cohort")
exclusion_right <- 0.95
exclusion_gap <- unit(5, "pt")
exclusion_line_offset <- unit(14, "mm")

grid.newpage()
flowchart(
  sources = list(
    ehr = boxGrob("Hospital EHR\nadmissions\nn = 241,820",
                  box_gp = source_gp),
    registry = boxGrob("Quality registry\nprocedures\nn = 38,420",
                       box_gp = source_gp),
    deaths = boxGrob("Population registry\nfollow-up\nn = 100%",
                     box_gp = source_gp)
  ),
  linked = boxGrob(
    "Linked study base\nunique patients with follow-up\nn = 29,614",
    x = main_x,
    box_gp = link_gp,
    width = unit(72, "mm")
  ),
  exclusions = list(
    prior = boxGrob("Previous diagnosis\nn = 4,108",
                    just = "left", box_gp = side_gp,
                    width = unit(42, "mm")),
    missing = boxGrob("Missing key\ncovariates\nn = 962",
                      just = "left", box_gp = side_gp,
                      width = unit(42, "mm")),
    outside = boxGrob("Outside study\nwindow\nn = 1,327",
                      just = "left", box_gp = side_gp,
                      width = unit(42, "mm"))
  ),
  cohort = boxGrob(
    "Primary cohort\nn = 23,217",
    box_gp = cohort_gp,
    width = unit(62, "mm")
  ),
  outputs = list(
    primary = boxGrob("Primary analysis\ncomplete case\nn = 22,144",
                      box_gp = final_gp),
    imputed = boxGrob("Sensitivity analysis\nmultiple imputation\nn = 23,217",
                      box_gp = final_gp),
    negative = boxGrob("Negative control\noutcome check\nn = 21,903",
                       box_gp = final_gp)
  )
) |>
  spread(axis = "y", margin = unit(8, "mm"), exclude = "exclusions") |>
  equalizeWidths(subelement = main_path) |>
  align(axis = "x", subelement = "cohort", reference = "linked") |>
  move(subelement = c("exclusions", "prior"),
       y = position("linked", position = "center", type = "y") - exclusion_gap,
       just = c(NA, "top")) |>
  move(subelement = c("exclusions", "missing"),
       y = position(c("exclusions", "prior"), position = "bottom", type = "y") - exclusion_gap,
       just = c(NA, "top")) |>
  move(subelement = c("exclusions", "outside"),
       y = position(c("exclusions", "missing"), position = "bottom", type = "y") - exclusion_gap,
       just = c(NA, "top")) |>
  equalizeWidths(subelement = "sources") |>
  equalizeWidths(subelement = "exclusions", width = unit(42, "mm")) |>
  equalizeWidths(subelement = "outputs") |>
  move(subelement = "exclusions", x = exclusion_right, just = "right") |>
  spread(axis = "x", subelement = "sources", margin = source_margin, type = "center") |>
  spread(axis = "x", subelement = "outputs", margin = output_margin, type = "center") |>
  connect("sources", "linked", type = "vertical_axis", lty_gp = con_gp, arrow_size = 3) |>
  connect("linked", "cohort", type = "v", lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  connect("linked", "exclusions",
          type = "side", lty_gp = excl_gp, arrow_size = 3,
          side = "right", end_side = "left",
          side_route = "outside",
          side_offset = exclusion_line_offset,
          label = "Excluded\nn = 6,397",
          label_gp = gpar(col = "#AD1457", cex = 0.8)) |>
  connect("cohort", "outputs", type = "N", lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  print()

The audit trail is part of
the story

Another common workflow is less glamorous but just as important: data
validation.

Many research projects have a small data-engineering pipeline even
when nobody calls it that. Data arrive through forms, imports, manual
entry, and collaborator spreadsheets. Then someone checks missing
fields, duplicates, impossible dates, inconsistent IDs, and
outliers.

That process is often hidden in prose. A compact flowchart can make
it visible without turning the methods section into a systems manual. It
is also a useful project-management figure: the same chart can be shown
to clinicians, data managers, statisticians, and co-authors.

Note how the box shapes carry meaning here — ellipses, databases,
documents, tapes, and diamonds all come from dedicated
box*Grob() helpers:

input_gp <- gpar(fill = "#F3F8FF", col = "#3B73C5", lwd = 1.3)
process_gp <- gpar(fill = "#FFF4C7", col = "#C69214", lwd = 1.3)
issue_gp <- gpar(fill = "#FCE4EC", col = "#AD1457", lwd = 1.2)
output_gp <- gpar(fill = "#E8F5E9", col = "#2E7D32", lwd = 1.3)
note_gp <- gpar(fill = "#FFFFFF", col = "#607D8B", lwd = 1, lty = 2)
con_gp <- gpar(col = "#555555", fill = "#555555", lwd = 1.3)
issue_con_gp <- gpar(col = "#AD1457", fill = "#AD1457", lwd = 1.1)

main_path <- list("validation", "clean")
issue_column_x <- 0.08
log_column_x <- 0.92
input_shape_width <- unit(42, "mm")
input_shape_height <- unit(24, "mm")
issue_shape_width <- unit(48, "mm")
issue_shape_height <- unit(14, "mm")

grid.newpage()
flowchart(
  inputs = list(
    web = boxEllipseGrob("REDCap\nform",
                         width = input_shape_width,
                         height = input_shape_height,
                         box_gp = input_gp),
    import = boxDatabaseGrob("CSV\nimport",
                             width = input_shape_width,
                             height = input_shape_height,
                             box_gp = input_gp),
    manual = boxDocumentGrob("Manual\nentry",
                             width = input_shape_width,
                             height = input_shape_height,
                             box_gp = input_gp)
  ),
  shape_note = boxGrob(
    "Shape indicates\nsource type",
    just = "left",
    width = unit(36, "mm"),
    box_gp = note_gp
  ),
  validation = boxTapeGrob(
    "Validation queue\nIDs, dates, ranges, missingness",
    width = unit(.58, "npc"),
    height = unit(.14, "npc"),
    box_gp = process_gp
  ),
  issues = list(
    missing = boxDiamondGrob("Missing\nfields",
                             width = issue_shape_width,
                             height = issue_shape_height,
                             box_gp = issue_gp),
    duplicate = boxDiamondGrob("Duplicate\nID",
                               width = issue_shape_width,
                               height = issue_shape_height,
                               box_gp = issue_gp),
    outlier = boxDiamondGrob("Outlier\nvalue",
                             width = issue_shape_width,
                             height = issue_shape_height,
                             box_gp = issue_gp)
  ),
  log = boxDocumentsGrob(
    "Issue log\nqueries sent\nchanges reviewed",
    width = unit(48, "mm"),
    height = unit(.44, "npc"),
    box_gp = issue_gp
  ),
  clean = boxDatabaseGrob(
    "Analysis-ready dataset\nlocked for report",
    width = unit(.44, "npc"),
    height = unit(.16, "npc"),
    box_gp = output_gp
  )
) |>
  spread(axis = "y", margin = unit(7, "mm"),
         exclude = list("issues", "shape_note")) |>
  spread(axis = "x", subelement = "inputs",
         from = 0, to = 0.7, margin = 0.05,
         type = "center") |>
  equalizeWidths(subelement = main_path) |>
  align(axis = "x", subelement = "validation", reference = "inputs") |>
  align(axis = "x", subelement = "clean", reference = "validation") |>
  align(axis = "y", subelement = "shape_note", reference = "inputs") |>
  align(axis = "y", subelement = "log",
        references = list("validation", "clean")) |>
  spread(axis = "y", subelement = "issues",
         from = position("log", position = "top", type = "y"),
         to = position("log", position = "bottom", type = "y"),
         margin = unit(2, "mm")) |>
  move(subelement = "shape_note", x = 0.95, just = "right") |>
  move(subelement = "issues", x = issue_column_x, just = "left") |>
  move(subelement = "log", x = log_column_x, just = "right") |>
  connect("inputs", "validation",
          type = "vertical_axis", lty_gp = con_gp, arrow_size = 3) |>
  connect("issues", "log",
          type = "horizontal_axis", lty_gp = issue_con_gp, arrow_size = 3) |>
  connect("validation", "clean", type = "vertical_axis",
          lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  print()

Follow-up is rarely just
down the page

Longitudinal studies often need to distinguish between people who are
lost, censored, withdrawn, dead, or still contributing information up to
a time point. A simple downward flow can imply that everyone leaving a
box disappears from the analysis, which is not always true.

Dotted return arrows are useful for this. They can show that a
participant left direct follow-up but still contributes information to
the final analysis up to censoring. That is a visual detail, but it
communicates an analytical idea.

This is where small flowchart improvements matter. Not because the
reader cares about the drawing API, but because the figure can express
the study design more faithfully.

options(boxGrobTxtPadding = unit(1, "mm"))

main_gp <- gpar(fill = "#FFFFFF", col = "#263238", lwd = 1.2)
arm_gp <- gpar(fill = "#E3F2FD", col = "#1565C0", lwd = 1.3)
ex_gp <- gpar(fill = "#FFF8E1", col = "#C69214", lwd = 1.2)
con_gp <- gpar(col = "#1565C0", fill = "#1565C0", lwd = 1.3)
side_gp <- gpar(col = "#C69214", fill = "#C69214", lwd = 1.2)
dotted_gp <- gpar(col = "#455A64", fill = "#455A64", lwd = 1.1, lty = 2)

arm_from <- .24
arm_to <- .76
box_width <- unit(54, "mm")
ex_width <- unit(45, "mm")
ex_page_margin <- 0.03           # excluded columns hug the page edge by this npc margin
side_offset <- unit(4, "mm")     # side branches step out this far before turning to the excluded box
fan_in_offset <- unit(2, "mm")   # dotted return line runs 2 mm outside the excluded boxes

grid.newpage()
flowchart(
  rando = boxGrob("Randomised\nN = 197", box_gp = main_gp),
  groups = list(
    boxGrob("96 assigned to intervention\n95 received treatment",
            box_gp = arm_gp),
    boxGrob("101 assigned to control\n93 received treatment",
            box_gp = arm_gp)
  ),
  ex1 = list(
    boxGrob("8 died\n1 withdrew consent", just = "left", box_gp = ex_gp),
    boxGrob("18 died\n1 withdrew consent", just = "left", box_gp = ex_gp)
  ),
  groups1 = list(
    boxGrob("87 completed day 30\nfollow-up", box_gp = arm_gp),
    boxGrob("79 completed day 30\nfollow-up", box_gp = arm_gp)
  ),
  ex2 = list(
    boxGrob("8 died", just = "left", box_gp = ex_gp),
    boxGrob("9 died\n1 withdrew consent\n2 lost to follow-up",
            just = "left", box_gp = ex_gp)
  ),
  groups2 = list(
    boxGrob("79 completed day 180\nfollow-up", box_gp = arm_gp),
    boxGrob("68 completed day 180\nfollow-up", box_gp = arm_gp)
  ),
  analysis = list(
    boxGrob("95 included in primary\noutcome analysis", box_gp = arm_gp),
    boxGrob("95 included in primary\noutcome analysis", box_gp = arm_gp)
  )
) |>
  spread(axis = "y", margin = unit(0.02, "npc")) |>
  equalizeWidths(subelement = stringr::regex("^groups|analysis"), width = box_width) |>
  equalizeHeights(subelement = stringr::regex("^groups|analysis")) |>
  equalizeWidths(subelement = stringr::regex("^ex"), width = ex_width) |>
  spread(subelement = stringr::regex("^groups|analysis"), axis = "x",
         from = arm_from, to = arm_to, type = "center") |>
  move(subelement = "rando",
       x = position("groups", position = "center", type = "x")) |>
  move(subelement = list(c("ex1", 1), c("ex2", 1)),
       x = ex_page_margin, just = "left") |>
  move(subelement = list(c("ex1", 2), c("ex2", 2)),
       x = 1 - ex_page_margin, just = "right") |>
  connect("rando", "groups", type = "N", lty_gp = con_gp, arrow_size = 3, smooth = TRUE) |>
  connect(c("groups$1", "groups1$1"), c("ex1$1", "ex2$1"),
          type = "side", lty_gp = side_gp, arrow_size = 3,
          side = "left", end_side = "right",
          side_route = "outside", side_offset = side_offset) |>
  connect(c("groups$2", "groups1$2"), c("ex1$2", "ex2$2"),
          type = "side", lty_gp = side_gp, arrow_size = 3,
          side = "right", end_side = "left",
          side_route = "outside", side_offset = side_offset) |>
  connect("groups", "groups1", type = "vertical", lty_gp = con_gp, arrow_size = 3) |>
  connect("groups1", "groups2", type = "vertical", lty_gp = con_gp, arrow_size = 3) |>
  connect("groups2", "analysis", type = "vertical", lty_gp = con_gp, arrow_size = 3) |>
  connect(list("ex1$1", "ex2$1"), "analysis$1", type = "side",
          lty_gp = dotted_gp, arrow_size = 3,
          side = "left", end_side = "left",
          side_route = "outside",
          side_offset = fan_in_offset) |>
  connect(list("ex1$2", "ex2$2"), "analysis$2", type = "side",
          lty_gp = dotted_gp, arrow_size = 3,
          side = "right", end_side = "right",
          side_route = "outside",
          side_offset = fan_in_offset) |>
  print()

Why this belongs in Gmisc

Gmisc has always collected the small tools I found myself needing
around medical statistics: descriptive tables, transition plots, and
grid-based figures. Flowcharts fit that same pattern. They are not a
statistical model, but they are part of how research is
communicated.

The new flowchart work in 3.4.0 is therefore aimed at the practical
problems:

• making CONSORT-like diagrams less painful to draw
• keeping grouped stages aligned and readable
• making arrows behave predictably
• supporting side paths, return paths, and repeated box patterns
• producing figures that can be regenerated when the study
  changes

The vignette contains the full API and examples:

vignette("Grid-based_flowcharts", package = "Gmisc")

The blog figures in this post are intentionally close to things
researchers already have in their workflow: trial enrollment, registry
construction, data validation, and follow-up accounting. My hope is that
they make the flowchart tools feel less like a drawing utility and more
like a small extension of the analysis itself.

EFAEFA explores unknown factor structures while CFA tests a predefined model. Compare both methods side by side with R code using psych and lavaan — know which to use in your dissertation.

Key Points

1. EFA and CFA are both factor analysis methods, but they serve opposite purposes: EFA discovers factor structure; CFA tests a pre-specified structure.
2. The core difference lies in factor loadings — EFA lets all items load freely on all factors; CFA constrains items to load only on their pre-assigned factor.
3. In R, EFA uses the psych package (fa()); CFA uses the lavaan package (cfa()).
4. CFA requires goodness-of-fit evaluation: CFI, TLI, RMSEA, and SRMR — all with established thresholds.
5. Many dissertations use both in sequence — EFA on a pilot sample, CFA on an independent main sample.

EFA vs CFA: Quick Comparison

EFA vs CFA: What Is the Difference Between Exploratory and Confirmatory Factor Analysis?

EFA and CFA are both forms of factor analysis — statistical methods that model the relationships between observed variables (e.g., survey items) and unobserved latent variables called factors. They are not competing methods; they serve different phases of measurement research. 

EFA (Exploratory Factor Analysis) lets the data reveal its own factor structure when you have no strong prior theory.

CFA (Confirmatory Factor Analysis) tests whether your observed data fit a structure you have already specified based on theory or prior EFA results.

Factor Anaysis

Dissertation quick-pick rule:
Using a new or adapted questionnaire with no established factor structure? → Start with EFA.
Replicating an established scale (Big Five, JSS, UTAUT) in a new sample? → Use CFA directly.
Both in the same study? → EFA on the pilot sample; CFA on the main independent sample.

What Is Exploratory Factor Analysis (EFA)?

Exploratory Factor Analysis (EFA) is a data-driven method that identifies the number and nature of latent factors underlying a set of observed variables, without imposing any prior constraints on which items load on which factors. EFA is theory-generating — it reveals patterns in your data that can later be formalised into a testable model for CFA.

Example: You design a 20-item questionnaire to measure academic motivation. You have no prior theory about how many sub-dimensions exist. EFA will cluster those 20 items into 3–5 factors (e.g., intrinsic motivation, extrinsic motivation, self-regulation) based purely on their intercorrelations — and the factor structure emerges from the data, not from your assumptions.

When to Use EFA

EFA Assumptions and Sample Size Requirements

Before running EFA, your data must meet several requirements. Items should be measured on interval or ordinal scales — Likert scales are appropriate. Check data factorability using the Kaiser-Meyer-Olkin (KMO) test (value > 0.60 required; > 0.80 is good) and Bartlett’s test of sphericity (p < 0.05 required). For sample size, the recommended minimum is 100 participants, but most methodologists advise at least 5–10 cases per item. A 20-item scale needs at least 200 respondents for stable EFA results.

How to Run EFA in R Using the psych Package

The psych package provides the most complete EFA workflow in R. The example below uses the built-in bfi dataset (25 Big Five personality items from the psych package itself):

# Step 1: Install and load the psych package
install.packages("psych")
library(psych)

# Step 2: Load data — bfi = Big Five Inventory (25 personality items)
data(bfi)
bfi_items <- bfi[, 1:25]          # Select only the 25 personality items

# Step 3: Test factorability before running EFA
KMO(bfi_items)                     # KMO > 0.60 required
cortest.bartlett(bfi_items, n = nrow(bfi_items))  # p < 0.05 required

# Step 4: Determine the number of factors via parallel analysis
fa.parallel(bfi_items, fm = "ml", fa = "fa")

# Step 5: Run EFA with 5 factors and oblimin (oblique) rotation
efa_model <- fa(bfi_items,
                nfactors = 5,
                rotate    = "oblimin",  # oblique — factors are allowed to correlate
                fm        = "ml")       # maximum likelihood estimation

# Step 6: Inspect factor loadings (show only loadings > 0.30)
print(efa_model$loadings, cutoff = 0.3)

# Step 7: View factor structure diagram
fa.diagram(efa_model)

Rotation choice:
Use rotate = "oblimin" (oblique) as your default — psychological and social science constructs are almost always correlated. Only use rotate = "varimax" (orthogonal) if you have a strong theoretical reason to assume completely independent factors.

What Is Confirmatory Factor Analysis (CFA)?

Confirmatory Factor Analysis (CFA) is a theory-testing method. You specify a measurement model in advance — exactly how many factors exist, which items load on which factors, and whether factors are correlated — then evaluate how well your observed data fit that model using goodness-of-fit indices. CFA is part of the Structural Equation Modelling (SEM) framework and is the standard method for establishing construct validity in dissertation research.

Example: Prior research and your literature review both support a 2-factor model of academic motivation: intrinsic and extrinsic motivation. You specify this 2-factor CFA model with 10 items assigned in advance, fit it to your data, and evaluate whether CFI > 0.95, RMSEA < 0.06, and SRMR < 0.08. If fit is acceptable, you have confirmed the structure and can proceed to SEM.

When to Use CFA

How to Run CFA in R Using the lavaan Package

The lavaan package is the standard CFA and SEM tool in R. Here is a full working example using two factors from the bfi dataset:

# Step 1: Install and load lavaan
install.packages("lavaan")
library(lavaan)

# Step 2: Define your measurement model using lavaan syntax
# Each line: factor_name =~ item1 + item2 + item3 ...
cfa_model <- '
  agreeableness     =~ A1 + A2 + A3 + A4 + A5
  conscientiousness =~ C1 + C2 + C3 + C4 + C5
'

# Step 3: Fit the model to your data
fit <- cfa(cfa_model,
            data   = bfi,
            std.lv = TRUE)   # standardise latent variables

# Step 4: Full model summary with standardised loadings and fit indices
summary(fit, fit.measures = TRUE, standardized = TRUE)

# Step 5: Extract specific fit indices for reporting
fitMeasures(fit, c("cfi", "tli", "rmsea", "rmsea.ci.lower",
                   "rmsea.ci.upper", "srmr", "chisq", "df", "pvalue"))

# Step 6: Inspect modification indices if fit is poor
modindices(fit, sort. = TRUE, maximum.number = 10)

EFA vs CFA: Full Head-to-Head Comparison

CFA Model Fit Indices: RMSEA, CFI, TLI, and SRMR Explained

When you run CFA, evaluating model fit is not optional — it is the core output. A CFA result without fit indices is unpublishable. The table below shows every index you need to report, what it measures, and the widely accepted thresholds. Report at least three; never rely on χ² alone (it is highly sensitive to sample size).

APA reporting template for dissertations:
"The two-factor CFA model demonstrated acceptable fit: χ²(34) = 67.2, p < .001, CFI = .96, TLI = .95, RMSEA = .047 [90% CI: .029–.064], SRMR = .051. All standardised factor loadings were statistically significant and exceeded .50 (range: .53–.78)."

Can You Use Both EFA and CFA in the Same Dissertation?

Yes — and in many quantitative dissertations, using both is the most rigorous approach. The critical rule: you must use independent datasets. Using the same data for EFA and then CFA is a methodological error that peer reviewers and examiners will flag, because a CFA model derived from EFA results will always fit the same data well — that is circularity, not validation.

The correct sequential approach, step by step:

1. Collect two independent datasets. Option A: run a pilot study (n ≥ 100–150) for EFA, then collect a main sample (n ≥ 200–300) for CFA. Option B: collect one large dataset and split it randomly 50/50.
2. Run EFA on Sample 1 using the psych package in R. Report KMO, Bartlett's test, parallel analysis output, factor loadings, communalities, and percentage variance explained.
3. Specify the CFA model based on the EFA factor structure. Assign each item to the factor it loaded most strongly on. Drop items with cross-loadings above 0.30 on two or more factors.
4. Run CFA on Sample 2 using lavaan. Fit the model, evaluate fit indices, and inspect modification indices if fit is inadequate.
5. Report both analyses in your methodology chapter, clearly stating which sample was used for which analysis and why the sequential approach was chosen.

How to Choose Between EFA and CFA: Decision Rules

The choice between EFA and CFA depends on your research question, the state of the literature, and the purpose of your factor analysis. These rules cover the most common dissertation scenarios:

• No prior theory about factor structure → EFA
• Established, well-cited factor structure from prior studies → CFA
• Adapting a scale to a new language, culture, or population → EFA first, then CFA
• Building toward Structural Equation Modelling → CFA mandatory
• Developing a new psychometric scale from scratch → EFA first, CFA to validate
• Comparing two competing theoretical models → CFA (use model comparison with likelihood ratio test)
• Mixed or contradictory literature on factor structure → EFA
• Testing construct validity (convergent + discriminant) → CFA

Examples of EFA and CFA in Research

Example 1: EFA of Personality Traits — The Big Five

The most influential application of EFA is the development of the Big Five personality model. Researchers applied EFA to hundreds of personality adjectives across multiple independent samples. Without any prior constraint on which traits should cluster together, EFA consistently revealed five factors: Openness, Conscientiousness, Extraversion, Agreeableness, and Neuroticism. This is EFA at its best — no prior theory constrained the analysis, and the five-factor structure replicated across cultures and languages.

To replicate this in R, run EFA with 5 factors and oblimin rotation on the bfi dataset in the psych package (code shown above). Each of the five resulting factors maps cleanly onto one of the Big Five dimensions, with factor loadings above 0.40 for the primary items.

Example 2: CFA of Job Satisfaction — The JSS

Spector's (1985) Job Satisfaction Survey (JSS) proposes a 9-factor model covering pay, promotion, supervision, fringe benefits, contingent rewards, operating procedures, co-workers, nature of work, and communication. A researcher validating the JSS in a healthcare sample would use CFA: specify all 36 items loading on their designated factors, fit the model with lavaan, and evaluate whether CFI > 0.95, RMSEA < 0.06, and SRMR < 0.08.

If fit is poor (e.g., RMSEA > 0.08), consult modification indices and consider whether any theoretically justifiable correlated residuals between items within the same facet would improve fit. Only modify parameters where there is both statistical and substantive justification. Need help interpreting your CFA output or writing up your p-values and fit indices? Message Dr. Zubair on WhatsApp.

Common EFA and CFA Mistakes to Avoid

• Using EFA and CFA on the same sample — the most common dissertation error. Always use independent samples.
• Choosing the number of EFA factors by eigenvalue > 1 rule alone — this systematically over-extracts factors. Use parallel analysis (fa.parallel()) instead.
• Using orthogonal rotation (varimax) by default — most psychological and social science constructs are correlated. Use oblimin rotation unless theory dictates independence.
• Reporting only χ² for CFA — χ² is significant with any n > 200. Always report CFI, TLI, RMSEA, and SRMR alongside it.
• Keeping cross-loading items in the CFA model — items that loaded > 0.30 on two or more factors in EFA should be dropped before CFA specification.
• Fewer than three items per factor — two-item factors are under-identified in CFA. Each factor needs at least 3 indicators.
• Confusing EFA with PCA — Principal Component Analysis (PCA) is a data reduction method, not a factor analysis technique. See our guide on Factor Analysis vs PCA for the full distinction.

EFA and CFA in Thesis and Dissertation: Reporting Requirements

Both analyses appear in the methodology chapter under "Measurement Validation" or "Scale Development." Here is exactly what your committee and journal reviewers expect:

• For EFA: Report KMO value, Bartlett's test (χ², df, p), number of factors extracted, method for determining factor number (parallel analysis recommended), extraction method (ML recommended), rotation method, eigenvalues for retained factors, percentage variance explained by each factor and total, and a complete factor loading matrix with items bolded if > 0.30.
• For CFA: Report the measurement model diagram (path diagram), sample size, estimation method (ML for continuous/normal data; WLSMV for ordinal/Likert), and fit indices: χ²(df), p-value, CFI, TLI, RMSEA [90% CI], SRMR. Report all standardised factor loadings and their significance. State whether modification indices were consulted and, if the model was modified, provide both theoretical and statistical justification.
• For both in the same study: Clearly label Sample 1 (EFA) and Sample 2 (CFA) in your methods section. Justify the sequential strategy. Provide descriptive statistics for both samples.

For related analyses your dissertation may also require, see our guides on EFA in R with psych, PCA in R, Factor Analysis vs PCA, and normality testing with Shapiro-Wilk before running your analyses.

Conclusion

The difference between EFA and CFA comes down to one question: do you know the factor structure, or are you trying to find it? EFA discovers structure from data when theory is absent or weak. CFA confirms a structure you have already specified when theory is strong or prior EFA results exist. In R, the psych package handles EFA and the lavaan package handles CFA — both are free, well-documented, and the current standard in academic research.

For PhD and Master's dissertation researchers: the most defensible methodology for a new measurement instrument is EFA on a pilot sample followed by CFA on an independent main sample. This sequential approach demonstrates both exploratory rigour and confirmatory validity to examiners and reviewers.

If you need expert help with your EFA or CFA analysis — including running the analysis in R or SPSS, interpreting fit indices, writing up results in APA format, or preparing your methodology chapter — contact Dr. Zubair Goraya on WhatsApp or book a session via the link below.

Frequently Asked Questions

With the recent resignation announcement from United Kingdom (UK) Prime Minister Keir Starmer, there have been a flurry of people talking about how many UK prime ministers there have been in the past decade, short terms for prime ministers, and so on. I wanted a historical perspective and so grabbed the data from Wikipedia. Wikipedia has a convenient single table list of all of the UK prime ministers since the term began being used informally by Robert Walpole. Walpole was effectively prime minister of the Kingdom of Great Britain from 1721 onwards.

The first prime minister of the United Kingdom of Great Britain and Ireland was William Pitt in 1801; and of United Kingdom of Great Britain and Northern Ireland was Andrew Bonar Law in 1922. But these distinctions will be largely disregarded for the purpose of this blog post.

Downloading prime ministerial data from Wikipedia

Here’s code to download and import that list from Wikipedia. This worked as at 23 June 2026, but Wikipedia pages are known to change in format so it’s brittle about whether it will work forever:

library(rvest)
library(tidyverse)
library(janitor)
library(slider) # for rolling sum
library(scales)
library(ggrepel)
library(kableExtra)

#-----------------Import and process data------------------------

url <- "https://en.wikipedia.org/wiki/List_of_prime_ministers_of_the_United_Kingdom"

page <- read_html(url)

# The main PM table is the first wikitable on the page
pm_table <- page |>
  html_element("table.wikitable") |>
  html_table(fill = TRUE) |> 
  clean_names() |> 
  select(
    pm = prime_minister_office_lifespan,
    start = term_of_office,
    end = term_of_office_2
  ) |> 
  # drop second line of column titles:
  slice(-1) |> 
  # find the PMs' names - everything up to the first [
  mutate(pm = str_extract(pm , ".*?\\["),
         pm = str_replace(pm, "\\[", ""),
        ) |> 
  # strip all the footnotes and stuff from the dates:
  mutate(across(everything(), ~ str_remove_all(.x, "\\[.*?\\]"))) |>  # remove [1] refs
  mutate(across(everything(), str_squish)) |> 
  mutate(start = as.Date(start, format = "%d %B %Y"),
         end = as.Date(end, format = "%d %B %Y"),
         end = if_else(is.na(end) & pm == "Keir Starmer",
                       as.Date("2026-07-10"),
                       end),
        duration = as.numeric(end - start),) |> 
  distinct() |> 
  mutate(pm = fct_reorder(pm, start, .desc = TRUE)) |> 
  group_by(pm) |> 
  mutate(last_end = max(end)) |> 
  ungroup()

Most and longest serving prime ministers

This lets us do some simple analysis. First, here are the UK prime ministers who have served the most often—that is, had more than one term:

Since the mid twentieth century, only Churchill and Wilson have had a second chance to be prime minister. In the nineteenth century it was much more common, with big names like Gladstone, Disraeli and Gascoyne-Cecil dominating politics while in government and out.

Here are those who have served the longest durations in total:

The first UK prime minister, Robert Walpole, was also the longest serving. From the 20th and 21st century, only Thatcher and Blair make the top ten list.

Those two simple tables were produced with this code:

#------------summary highlights----------

# prime ministers number of terms and total duration:
pm_summary <- pm_table |> 
  rename(`Prime minister` = pm) |> 
  group_by(`Prime minister`) |> 
  summarise(Terms = length(`Prime minister`),
            `Earliest start` = min(start),
            `Latest finish` = max(end),
            `Total duration` = sum(duration)) |> 
  arrange(desc(Terms), `Earliest start`) 

# Prime ministers with more than one term:
pm_summary |> 
  filter(Terms > 1) |> 
  kable() |> 
  kable_styling() 

# Longest serving prime ministers:
pm_summary |> 
  arrange(desc(`Total duration`)) |> 
  slice(1:10) |> 
  kable() |> 
  kable_styling() 

Graphic summaries

Tables are nice but graphics are better. Here is my attempt to summarise all the prime ministers of the UK (and of the predecessor Kingdom of Great Britain) in one picture. You probably need a full-sized screen for this, but with the right display I think the Gantt chartish style works nicely.

That chart produced with this code. There are a few clutter-minimisation polishing details here on top of my usual blog style, like suppressing the y axis labels and adding them instead as text close to the data. Much easier to read. And suppressing the horizontal gridlines.

--------------Draw plots-------------
the_title <- "Prime ministers of the United Kingdom and its predecessors, 1721 to 2026"

pm_table |> 
  ggplot(aes(y = pm, yend = pm)) +
  geom_segment(aes(x = start, xend = end),
               linewidth = 2, colour = "steelblue") +
  geom_text(data = distinct(pm_table, pm, last_end),
            aes(label = pm, x = last_end + 500),
            size = 2, hjust = 0, colour = "grey50") +
  scale_x_date(
    breaks = seq(as.Date("1720-01-01"), as.Date("2035-01-01"), by = "20 years"),
    date_labels = "%Y",
    sec.axis = sec_axis(~.),
  ) + 
  labs(x = "Year",
       y = "",
       title = the_title) +
  theme(axis.text.y = element_blank(),
        panel.grid.major.y  = element_blank(),
        panel.border = element_blank(), 
        axis.ticks.y = element_blank())

Secondly, it seems highly relevant to produce a plot of the distribution of durations:

And one showing the trend (or lack of trend) in durations over time:

Those two simple plots produced with this code. Perhaps the only point of particular interest here is how I used a subset of the data to highlight the names of prime ministers with term duration of less than 120 days or more than 3,000:

pm_table |> 
  ggplot(aes(x = duration)) +
  geom_density(colour = "steelblue") +
  geom_rug(colour = "steelblue") +
  scale_x_continuous(label = comma) +
  labs(x = "Duration in days",
       title = the_title)

pm_table |> 
  ggplot(aes(x = start, y = duration)) +
  geom_smooth(method = "gam", colour = "white") +
  geom_point(colour = "steelblue") +
  geom_text_repel(data = filter(pm_table, duration < 120 | duration > 3000), 
                   aes(label = pm), size = 2.8, seed = 123) +
  scale_y_sqrt(breaks = c(0.5, 1, 1:4 * 2) * 1000, label = comma) +
  labs(x = "Starting date of premiership",
       y = "Duration in days",
       title = the_title,
      subtitle = "Durations shown are of individual periods in office, not lifetime totals.")

Finally, the big question that seems to get a lot of attention. How many prime ministers per decade? Below is my effort at calculating and presenting this.

We can see that we are indeed going through a decade that is rich in UK prime ministers (and will be richer still in a month or so). But it’s not unprecedented. We’ve been at similar levels seeral times in the past, and in the politically turbulent 1830s there were even more premierships.

In fact, in the late twentieth century with Thatcher and Blair, the UK faced a period of unusually slow turnover of prime ministers. But that was a formative period in the life of many of today’s political commentators, so its not surprising that the current rapid turnover comes across as a surprise.

Code for this is below. Note that I calculated this on a daily basis. I’m not 100% I’ve got it right, but it passes my simplest reality checks (eg manually counting those we’ve had in the past ten years – six so far, although expeted soon to become seven).

cumulative_pms <- pm_table |> 
  full_join(tibble(start = seq(from = min(pm_table$start), 
                               to = max(pm_table$end), 
                               by = "1 day"))) |> 
  arrange(start) |> 
  mutate(starting_pms = if_else(is.na(pm), 0 , 1),
         rolling_pms = slide_sum(starting_pms, before = 3653),
         # I'm not sure I've got this right yet, but the idea is that in any given day,
         # the number of PMs in thepast 10 years is however many started in those 10 years,
         # plus 1 PM that you came into the period with. The exception being the time
         # of the very first prime minister, for which time you only have the rolling sum
         # of PMs that started:
         rolling_pms = if_else(start < (pm_table[1, ]$start + 3654), 
                               rolling_pms, 
                               rolling_pms + 1))

# When was the peak number of PMs in the last decade:
arrange(cumulative_pms, desc(rolling_pms))


cumulative_pms |> 
  ggplot(aes(x = start, y = rolling_pms)) +
  geom_line(colour = "steelblue") +
  scale_y_continuous(breaks = 0:max(cumulative_pms$rolling_pms)) +
  labs(x = "",
       y = "Number of prime ministers in past 10 years",
       title = the_title,
       subtitle = "Peak prime ministers per decade was in the 1830s")

That’s all for now.

How to Perform PCA in R: A Step-by-Step Tutorial Using prcomp()

To perform PCA in R, use the built-in prcomp() function: pca <- prcomp(data, scale = TRUE), then run summary(pca) to see how much variance each component explains. prcomp() handles centering, scaling, and the underlying singular value decomposition for you, and returns the loadings (pca$rotation) and component scores (pca$x) you need for any further analysis or visualization.

Complete Code — Copy, Paste, Run

# Complete PCA in R script -- copy, paste, run top to bottom

# 1. Install and load required packages
install.packages(c("ggplot2", "factoextra"))
library(ggplot2)
library(factoextra)
#2. Load the data and check assumptions
data(USArrests)
str(USArrests)
cor(USArrests)
apply(USArrests, 2, shapiro.test)
#3. Scale the data (prcomp() can also do this inline with scale = TRUE)
data_scaled <- scale(USArrests)
head(data_scaled)
#4. Run PCA
pca <- prcomp(USArrests, scale = TRUE)
summary(pca)
pca$rotation
#5. Decide how many components to keep
plot(pca$sdev^2, type = "b", xlab = "Principal Component", ylab = "Eigenvalue", main = "Scree Plot")
fviz_eig(pca, addlabels = TRUE, barfill = "steelblue", barcolor = "steelblue",
         linecolor = "firebrick", main = "Scree Plot: Variance Explained by Component")
pve <- pca$sdev^2 / sum(pca$sdev^2)
cum_pve <- cumsum(pve)
plot(cum_pve, type = "b", xlab = "Principal Component",
     ylab = "Cumulative Proportion of Variance Explained")
#6. Visualize results
biplot(pca, scale = 0)
fviz_pca_biplot(pca,
                geom.ind = "point",
                col.ind = as.factor(state.region),
                palette = "jco",
                addEllipses = TRUE,
                col.var = "black",
                repel = TRUE,
                legend.title = "Region") +
  theme_minimal()

Key Takeaways

• Principal component analysis (PCA) in R turns a set of correlated variables into a smaller set of uncorrelated principal components that capture most of the original variance.
• The built-in prcomp() function is the preferred way to run PCA in R because it uses singular value decomposition (SVD), which is numerically more stable than the older princomp() function.
• On the classic USArrests dataset, the first two principal components explain 86.75% of the total variance — verified directly from R output further down this page.
• You can decide how many components to keep using a scree plot, the Kaiser rule (eigenvalue > 1), or a cumulative variance threshold — this tutorial shows all three and where they disagree.
• PCA is not the same as Factor Analysis. The comparison table below explains exactly when to use each one.

If you have a dataset with many correlated variables and want to find the ones that actually matter, reduce its complexity, and visualize it in two dimensions without losing the patterns hidden inside it — PCA in R is the tool for the job.

My name is Dr Zubair Goraya. I hold a PhD-level background in statistics and have used R for statistical consulting and research for several years. I ran into the same questions you probably have now — how many components to keep, how to read the prcomp() output, how to build a biplot with ggplot2 — while working through PCA for my own research, and I have since helped thesis and dissertation students work through the exact same problem.

What Is Principal Component Analysis (PCA)?

Principal component analysis is a statistical technique for dimensionality reduction. It takes a set of variables that are correlated with one another and re-expresses them as a smaller set of new, uncorrelated variables called principal components. The first principal component (PC1) captures the largest possible amount of variance in the data; the second (PC2) captures the largest amount of the variance left over after PC1, and so on, with each component orthogonal to the ones before it.

How PCA Works Under the Hood

You don't need to compute this by hand — prcomp() does it for you — but understanding the four steps helps you interpret the output correctly:

1. Standardize each variable (mean 0, standard deviation 1) so no single variable dominates because of its scale.
2. Compute the covariance matrix of the standardized variables, which captures how every pair of variables moves together.
3. Decompose that matrix into eigenvectors (the direction of each principal component) and eigenvalues (how much variance that direction explains). The eigenvector with the largest eigenvalue is PC1.
4. Project the original data onto the new eigenvector axes. This produces the component scores — the new coordinates for every observation in the reduced space.

prcomp() performs this using singular value decomposition (SVD) rather than direct eigendecomposition of the covariance matrix, which is the same result reached through a numerically more stable route — this is why prcomp() is preferred over princomp() in R.

What You Need Before Running PCA in R

PCA only works on numeric, continuous variables. You'll use:

• prcomp() — built into base R's stats package, no installation needed.
• ggplot2 — for the custom scree plot and biplot in this tutorial.
• factoextra — a ggplot2-based wrapper that turns PCA output into publication-ready plots with a single function call.

install.packages(c("ggplot2", "factoextra"))
library(ggplot2)
library(factoextra)

This tutorial uses USArrests, a dataset built into R covering violent crime rates and urban population percentage for all 50 US states in 1973 — the same dataset used in An Introduction to Statistical Learning, so you can cross-check your output against a well-known reference.

Step 1: Check Your PCA Assumptions in R

PCA assumes your variables are numeric, continuous, linearly related, and reasonably normally distributed. Confirm this before you run anything:

data(USArrests)
str(USArrests)                       # numeric and continuous?
cor(USArrests)                       # linearly related?
apply(USArrests, 2, shapiro.test)    # normally distributed?

Running a Shapiro-Wilk normality test on each variable shows that UrbanPop is the only variable that does not significantly depart from normality (p = 0.977); Murder, Assault, and Rape all return p-values below 0.05:

PCA is fairly robust to mild non-normality, especially with n = 50, so this tutorial proceeds with all four variables rather than dropping any. If your own data fails normality more severely, consider a transformation first, or revisit whether PCA is the right tool at all — see the PCA vs. Factor Analysis comparison further down this page.

Step 2: Scale and Center Your Data

Murder, Assault, and Rape are measured per 100,000 residents while UrbanPop is a percentage — different scales entirely. Without scaling, Assault (values up to 337) would dominate the first component purely because of its larger numbers, not because it is more important. prcomp() has a built-in scale. argument that handles this for you, but you can also normalize the data manually first using Z-score standardization if you want to inspect the scaled values:

data_scaled <- scale(USArrests)
head(data_scaled)

Step 3: Run PCA in R With prcomp()

Now run PCA on the full four-variable dataset, with scaling handled inline:

pca <- prcomp(USArrests, scale = TRUE)
summary(pca)

This is the actual summary(pca) output for USArrests:

PC1 alone explains 62.0% of the total variance; PC1 and PC2 together explain 86.75%. That means you can reduce four variables down to two principal components and still retain most of the information in the original data.

Looking at pca$rotation shows what each component represents: PC1 loads heavily and negatively on Murder (-0.536), Assault (-0.583), and Rape (-0.543), with a smaller loading on UrbanPop (-0.278) — so PC1 is best read as an overall violent crime axis. PC2 loads almost entirely on UrbanPop (0.873), making it an urbanization axis largely independent of crime rate.

Step 4: Decide How Many Components to Keep

There is no single universal rule. In practice, researchers triangulate using three methods — and on this dataset, they don't all agree, which is itself a useful teaching example.

Method 1: The Scree Plot

plot(pca$sdev^2, type = "b", xlab = "Principal Component", ylab = "Eigenvalue", main = "Scree Plot")

Look for the "elbow" — the point where the line flattens out. The components before the elbow are worth keeping; the ones after add little.

The ggplot2 / factoextra Version

For a publication-ready scree plot built on ggplot2 instead of base R graphics, factoextra gives you a single-line solution:

fviz_eig(pca, addlabels = TRUE, barfill = "steelblue", barcolor = "steelblue",
         linecolor = "firebrick", main = "Scree Plot: Variance Explained by Component")

fviz_eig() returns a standard ggplot object, so you can layer on any further ggplot2 theming (+ theme_minimal(), custom labels, color palettes) exactly as you would with any other ggplot2 chart.

Method 2: The Kaiser Rule (Eigenvalue > 1)

The Kaiser criterion says: keep any component whose eigenvalue (standard deviation squared) exceeds 1. On this dataset:

This is a good example of why the Kaiser rule shouldn't be used in isolation: PC2's eigenvalue (0.990) is essentially 1, so a strict cutoff would discard a component that still explains nearly a quarter of total variance and has a clear, interpretable meaning (urbanization). Most applied researchers would retain PC2 here despite the Kaiser rule's technical "no."

Method 3: Cumulative Variance Threshold

pve <- pca$sdev^2 / sum(pca$sdev^2)
cum_pve <- cumsum(pve)
plot(cum_pve, type = "b", xlab = "Principal Component",
     ylab = "Cumulative Proportion of Variance Explained")

A common threshold is 80–90% cumulative variance. Here, two components clear 80% (86.75%), so retaining PC1 and PC2 is the most defensible choice across all three methods combined.

Step 5: Visualize PCA Results — Biplot in R

A biplot overlays the component scores (observations) and the loadings (original variables) on the same PC1/PC2 plane.

Base R Biplot

biplot(pca, scale = 0)

PCA Biplot With ggplot2 (factoextra)

The base R biplot is functional but not very customizable. For a ggplot2 biplot with colored points, confidence ellipses, and repelled labels, use fviz_pca_biplot():

fviz_pca_biplot(pca,
                geom.ind = "point",
                col.ind = as.factor(state.region),
                palette = "jco",
                addEllipses = TRUE,
                col.var = "black",
                repel = TRUE,
                legend.title = "Region") +
  theme_minimal()

Because fviz_pca_biplot() returns a ggplot object, every ggplot2 layer works: swap palette, add + labs(title = "..."), or change the theme without touching the underlying PCA math.

Southern states cluster toward high values on PC1 (higher violent crime), while Western states spread further along PC2 (higher urbanization) — exactly the two axes the loadings predicted in Step 3.

PCA vs. Factor Analysis: What's the Difference?

Students frequently confuse PCA with Factor Analysis because both reduce a large set of variables down to a few. They are not interchangeable, and using the wrong one is a common mistake in thesis methodology chapters:

If you're building a scale to measure a psychological construct (anxiety, job satisfaction, brand loyalty), you almost certainly want Exploratory Factor Analysis, not PCA — even though both are run on similarly structured survey data.

How to Report PCA Results in APA Style

For a thesis or journal manuscript, your methods and results sections should report PCA in a standard format. Based on the USArrests output above, here is a template you can adapt:

Sample write-up:
"A principal component analysis (PCA) was conducted on four crime-related variables (Murder, Assault, UrbanPop, Rape) using the prcomp() function in R (R Core Team). Variables were standardized prior to analysis. Using the Kaiser criterion (eigenvalue > 1) in combination with a scree plot and an 80% cumulative variance threshold, two components were retained. PC1 explained 62.0% of the variance and was driven primarily by Murder, Assault, and Rape (loadings of -0.54, -0.58, and -0.54, respectively), suggesting a general violent crime dimension. PC2 explained an additional 24.7% of the variance and loaded almost exclusively on UrbanPop (0.87), representing an urbanization dimension. Together, the two components accounted for 86.75% of total variance."

Always include a loadings table (as shown in Step 3) and either a scree plot or a cumulative-variance plot as a figure when reporting PCA in a formal write-up — reviewers and committee members expect to see the basis for your retention decision, not just the final component count.

Further Analysis: Using PCA Output Downstream

Once you have component scores (pca$x), they become inputs for other analyses — for example, as features for k-means clustering on a reduced, decorrelated dataset, or to detect multicollinearity before fitting a regression model.

Conclusion

You've now run a complete PCA in R: checking assumptions, scaling the data, fitting prcomp(), choosing how many components to keep using three different methods, visualizing results in both base R and ggplot2, distinguishing PCA from Factor Analysis, and reporting the output in APA format. The same workflow applies directly to your own dataset — just swap USArrests for your data frame.

Frequently Asked Questions

For most users, the built-in prcomp() function (stats package) is sufficient and is preferred over princomp() for its numerical stability. For richer visualization, pair it with factoextra. For multivariate exploratory work beyond PCA — including correspondence analysis — FactoMineR is the more comprehensive package.

Use plot(pca) for a quick scree plot, biplot(pca) for a base R biplot of scores and loadings, or factoextra::fviz_pca_biplot(pca) for a fully customizable ggplot2 version with color, ellipses, and labels, as shown in Step 5 above.

prcomp() returns an object of class "prcomp" containing: sdev (standard deviations of each component), rotation (the loadings matrix), x (the component scores for each observation), and center/scale (the centering and scaling values used).

princomp() uses spectral (eigen) decomposition of the covariance matrix. prcomp() uses singular value decomposition (SVD) of the data matrix directly, which has slightly better numerical accuracy. R's own documentation recommends prcomp() for this reason.

There's no universal rule. Triangulate using a scree plot (look for the elbow), the Kaiser criterion (eigenvalue > 1), and a cumulative variance threshold (commonly 80–90%) — Step 4 above shows all three applied to the same dataset, including a case where they disagree.

Standard PCA requires numeric, continuous variables. For categorical data, use Multiple Correspondence Analysis (MCA) instead, available via FactoMineR::MCA(), which applies the same dimensionality-reduction logic to a contingency table of category frequencies.

The three most frequent errors: (1) skipping scale = TRUE, which lets large-magnitude variables dominate the first component; (2) not checking for missing values or outliers beforehand, which distorts the covariance matrix; and (3) misreading the scree plot elbow or ignoring loading signs when interpreting what a component represents.

PCA is widely used for image compression, facial recognition feature extraction, gene expression analysis in bioinformatics, customer segmentation in marketing, and as a preprocessing step before clustering or regression to remove multicollinearity.

If you're working through PCA, Factor Analysis, or any other multivariate technique for a thesis or dissertation and want a second pair of eyes on your output, get in touch below.

May was Open Source Software Maintainer Month. Behind every R package there is at least one person who responds to issues, reviews pull requests, keeps up with dependency changes, and makes sure everything still works. During Maintainer Month we wanted to celebrate rOpenSci’s package maintainer community.

The social media campaign

One of our commitments to our community is to amplify the people who make it work. Social media is one of the ways we do that, so we thought Maintainer Month would be a great opportunity to highlight the people behind the packages through a social media campaign.

To run this campaign, we first needed permission from our maintainers to feature them. In our annual maintainer survey, we asked whether they would be interested in being featured in a public spotlight, and many said yes.

We also reached out to current and past Champions from our Champions Program, which trains and supports R developers from historically underrepresented groups in the open science community.

The result was a month-long series of spotlights: one maintainer at a time, each card sharing who they are, where they come from, and what they maintain.

This campaign brought together 37 maintainers from 15 countries, maintaining more than 50 packages that together serve thousands of researchers and data practitioners around the world.

The diversity of this group reflects the diversity of the rOpenSci community: archaeologists, bioinformaticians, ecologists, economists, statisticians, sociologists, professors, PhD students, engineers and educators.

We created 39 posts on our accounts on LinkedIn and Mastodon, which is bridge to BlueSky. All the posts were shared by other people and organizations and received comments from grateful users.

Meet all 37 maintainers

Here is the full list of maintainers we celebrated in May.

• Alex Koiter maintains {mbquartR}, for working with Manitoba’s quarter-section land survey system in watershed and land management research.

• Andrea Gomez Vargas maintains {ARcenso}, for accessing and analyzing Argentina’s national census data in R. Champions project.

• Austin Koontz maintains {SymbiotaR2}, an R interface to the Symbiota platform for accessing and managing biodiversity occurrence data from natural history collections.

• Bilikisu Wunmi Olatunji maintains {chartkickR}, an R wrapper for the Chartkick JavaScript library that makes it easy to create beautiful interactive charts and visualizations from R. Champions project.

• Carolina Pradier maintains {eph}, for downloading and analyzing microdata from Argentina’s Permanent Household Survey, supporting labour and socioeconomic research. Champions project.

• Daniel Vartanian maintains {mctq}, for processing data from the Munich Chronotype Questionnaire in sleep and chronobiology research.

• Erick Navarro Delgado maintains {RAMEN}, for identifying associations between environmental exposures and molecular outcomes in multi-omics research. Champions project.

• Erika Siregar maintains {rplaywright}, an R interface to Microsoft Playwright for browser automation and web testing. Champions project.

• Ezekiel Adebayo Ogundepo maintains {bulkreadr}, for simplifying the bulk import of multiple files into R across a range of formats. Champions project.

• Francesca Palmeira maintains {pcir}, for modeling species interaction data and food web structures in conservation research. Champions project.

• Guadalupe Pascal maintains {matildaNLP}, a package with a specialized corpus of Spanish texts from the Matilda initiative to support research on gender-aware language processing and policy. Champions project.

• Haydée Svab maintains {odbr}, for accessing open data urban mobility from a series of cities in Brazil. Champions project.

• Jeroen Ooms maintains {magick}, {pdftools}, and {gert}, packages for image processing, PDF manipulation, and Git operations in R.

• Jonathan Keane maintains {dittodb}, which makes testing database-backed code easy by recording and replaying real database interactions so tests can run without a live connection.

• Julia Silge maintains {qualtRics}, for importing survey data from the Qualtrics platform directly into R.

• Karl Broman maintains {chromer} and {aRxiv}, for accessing chromosome data and the arXiv preprint server.

• Maëlle Salmon maintains {saperlipopette}, {babelquarto}, and {babeldown}, tools for learn how to use git, create multilingual Quarto documents, and support translations workflows.

• Marcelo S. Perlin maintains {yfR}, for importing financial data from Yahoo Finance into R.

• Marcos Prunello maintains {karel}, a package that brings the Karel the Robot programming environment to R, designed to teach programming concepts and computational thinking to beginners. Champions project.

• Mark Padgham maintains {pkgcheck}, which automates software checks for packages submitted to rOpenSci peer review.

• Mauro Loprete maintains {metasurvey}, for processing and analyzing household survey microdata using a metadata-driven approach. Champions project.

• Micha Silver maintains {rOPTRAM}, implementing the OPtical TRApezoid Model for estimating soil moisture from satellite imagery.

• Moritz Hennicke maintains {nuts}, for working with the EU’s Nomenclature of Territorial Units for Statistics, useful in regional economics and policy research.

• Pao Corrales maintains {agroclimatico}, for calculating agroclimatic indices and bioclimatic variables for agricultural and environmental research. Champions project.

• Peter Desmet maintains {frictionless}, for working with open data standards and publishing datasets.

• Philippe Massicotte maintains {rnaturalearth}, {rnaturalearthdata}, and {gitignore}, for working with natural earth map data and project utilities.

• Sam Albers maintains {tidyhydat}, for accessing Canadian hydrometric data in a tidy format.

• Steffi Lazerte maintains {weathercan}, for downloading Canadian weather data directly from Environment and Climate Change Canada.

• Tad Dallas maintains {helminthR}, for accessing the London Natural History Museum’s host-parasite database.

• Ronald M. Visser maintains {dendroNetwork}, for creating and analyzing networks in dendrochronological research, combining archaeology and data science.

• Sehrish Kanwal maintains {RNAsum}, for summarising and visualising RNA-seq data analysis results in clinical cancer genomics workflows. Champions project.

• Victor Ordu maintains {naijR}, a package of tools and utilities for working with data and maps about Nigeria. Champions project.

• Will Gearty maintains {rredlist}, for accessing IUCN Red List data on threatened species.

• Will Landau maintains {targets}, a pipeline toolkit that makes data analysis in R faster and fully reproducible by tracking dependencies and only re-running what has changed.

• Will Pearse maintains {suppdata}, for downloading supplementary data files directly from published scientific articles across major journals.

• Yi-Chin Sunny Tseng maintains {bbsTaiwan}, for accessing and analyzing data from Taiwan’s Breeding Bird Survey. Champions project.

• Zhian Kamvar maintains {tinkr}, for reading and writing Markdown documents in R as XML.

Thank you Maintainers!

Maintaining open source software is an act of generosity. It takes time that could be spent elsewhere, and it often goes unacknowledged. Every bug fix, every answered issue, every new feature and update is a small gift to the people who depend on that package.

We are grateful to all the rOpenSci maintainers. If you use any of these packages, consider saying thank you. You can also let us know how you use these packages by sharing your use case, that we will feature in our website.

Want to learn more? Explore the rOpenSci’s packages in our website and check all the other packages universes in R-Universe.

The European Bioconductor Conference 2026 (EuroBioC2026) took place from June 3-5, 2026, in Turku, Finland. Hosted by the University of Turku and the Finnish Society for Bioinformatics at BioCity, the conference brought together the Bioconductor community to showcase the latest developments in Bioconductor software packages and discuss emerging technologies shaping computational biology. This year’s conference welcomed 147 in-person participants from 23 countries. Across three days, attendees participated in keynote lectures, short and flash talks, workshops, poster sessions, Birds-of-a-Feather discussions, and community events. The conference also marked an important milestone for the project as Bioconductor celebrated its 25th anniversary. The figures below summarise EuroBioC2026 at a glance: 147 attendees from 23 countries, 4 keynote speakers, 25 speakers, 68 posters, 6 workshops, 9 flash talks, and 3 Birds-of-a-Feather sessions.

Preconference

Ahead of the main conference, EuroBioC2026 hosted two preconference events on June 1-2. These were delivered in collaboration with the University of Turku, CompLifeSci, the Finnish Society for Bioinformatics, and members of the Bioconductor community. Running in parallel over two days, the events allowed participants to either strengthen their analytical skills through hands-on training or contribute directly to the development of Bioconductor software through collaborative coding projects.

Workshop: Orchestrating Microbiome Analysis with Bioconductor

The preconference workshop focused on microbiome data analysis using Bioconductor and followed the Bioconductor Carpentry model, combining interactive instruction with practical exercises. Over two days, participants learned how to import, process, and analyse microbiome datasets using the established Bioconductor workflow; Orchestrating Microbiome Analysis (OMA). The workshop covered diversity analyses, differential abundance testing, and approaches for integrating microbiome data with other omics data types. A key component of the workshop was the use of cloud computing resources from CSC (Finnish IT Centre for Science) using Noppe, which provided participants with immediate access to all required datasets, software, and computing resources. By removing installation and configuration barriers, instructors were able to begin teaching immediately and spend more time focusing on the workshop content rather than troubleshooting technical issues. The platform also ensured that all participants worked within the same environment, creating a smoother learning experience for everyone involved. The workshop was hands-on throughout: participants asked questions, worked through exercises, and discussed how the workflows related to their own projects. This practical, open, instructor-led format aligns well with the approach shared by both Bioconductor and The Carpentries. The workshop was led by Leo Lahti, Himel Mallick, Thomaz Bastiaanssen, Tuomas Borman, and Giulio Benedetti.

Participants during the pre-conference microbiome workshop.

Hackathon

We held the first of our series of hackathons attached to Bioconductor conferences this June at EuroBioC2026 in Turku, Finland. Eighteen in-person attendees worked on four projects focused on interoperability, and at least three of those are now being prepared for submission to BioHackrXiv.

A big congratulations to all the participants for their effort. You can read more about the projects from the EuroBioC2026 Hackathon. We’re looking forward to building on this work at the North American Bioconductor conference, BioC2026, in Seattle this August. See the BioC2026 Hackathon for more details.

Participants during the pre-conference hackathon.

Programme overview

EuroBioC2026 covered both established and emerging areas of computational biology, with a consistent focus on reproducible and open-source research.

Keynotes

Keynotes at EuroBioC2026 covered functional genomics, machine learning, microbiome research, and the direction of computational biology. Across four talks, speakers addressed how data science is changing biological research, and what that means for reproducibility, interpretation, and open-source software.

Helena Kilpinen opened the conference with Morphological profiling of in vitro neurons: Visualizing complexity in cellular disease models. Her talk explored how high-content imaging and morphological profiling can be used to better understand cellular phenotypes in disease models. By combining large-scale imaging data with computational approaches, she demonstrated how researchers can uncover subtle cellular differences that may provide insights into disease mechanisms.

Anders Krogh presented A Deep Generative Model for Gene Expression and Multimodal Data, showcasing how modern machine learning approaches can be used to model increasingly complex biological datasets. His keynote highlighted the potential of generative models to integrate multiple data modalities and improve our understanding of gene regulation and cellular states.

Aura Raulo delivered a keynote titled Modeling the spread of microbial communities in contact networks. Drawing on concepts from ecology, microbiology, and network science, they explored how microbial communities are transmitted between individuals and populations. The talk examined how host-associated microbiomes are shaped and shared, and what drives their spread.

The final keynote was delivered by Levi Waldron, who addressed a topic now central to many scientific discussions: Bioconductor in the age of AI. What do we do now? His talk examined the opportunities and challenges that AI presents for open-source scientific software. He encouraged the community to think about how AI tools can complement existing work, without compromising the transparency, reproducibility, and scientific rigour the project has built over 25 years.

Short talks and flash talks

The short talks at EuroBioC2026 reflected the diversity of the Bioconductor community, spanning topics from single-cell and spatial biology to microbiome research, proteomics, metabolomics, and multi-omics data integration. Several presentations introduced new software packages and statistical methods aimed at improving reproducibility, scalability, and interoperability in biological data analysis. Alongside methodological advances, speakers also covered broader community topics, including sustainable open-source software, environmentally conscious computing, training initiatives, and the growing role of artificial intelligence in computational biology. Together, the talks provided a good picture of the scientific questions being addressed with Bioconductor and the people driving its development.

Poster sessions

The 68 posters presented at EuroBioC2026 covered a broad mix of biological applications and software development. Topics included spatial omics, microbiome research, proteomics, metabolomics, disease modelling, and machine learning, as well as new packages and infrastructure projects from across the Bioconductor ecosystem. The poster sessions encouraged interactions between package developers, researchers, students, and first-time conference attendees, helping strengthen collaborations across the community.

EuroBioC2026 participants during a poster session.

Birds-of-a-Feather sessions

The three 90-minute Birds-of-a-Feather (BoF) sessions offered attendees an opportunity to connect around shared interests and exchange experiences, discuss challenges, and share ideas. The sessions were proposed by participants during the conference including sessions focused on strengthening the Finnish Bioconductor community, supporting early-career researchers, and embedding environmental sustainability into Bioconductor packages and research workflows. One outcome from the early-career researcher discussion was the creation of a dedicated student–ECR Zulip channel to support continued connection within the community. The BoF sessions continued a tradition of community-led discussion that has been part of Bioconductor events for years.

Workshops

The workshop sessions offered attendees an opportunity to explore a range of Bioconductor tools and workflows through hands-on demonstrations led by community members. Topics included proteomics data analysis, integrative analysis of histopathological images and multi-omics data, ChIP-seq analysis, differential expression analysis, post-translational modification analysis, and interoperable mass spectrometry workflows combining R and Python. Participants had the opportunity to engage directly with instructors, ask questions, and learn how the presented tools could be applied to their own research projects. Together, the workshops showcased the breadth of analytical domains supported by the Bioconductor ecosystem.

Celebrating 25 years of Bioconductor

A major highlight of EuroBioC2026 was the celebration of Bioconductor’s 25th anniversary. Since its founding in 2001, Bioconductor has grown from a small collection of software packages into a global open-source community used by thousands of researchers worldwide. Over the past quarter-century, it has become a central resource for reproducible computational biology, providing infrastructure, software, training, and community support across numerous biological disciplines.

One of the highlights of the celebration was a retrospective presented by Maria Doyle, Bioconductor Community Manager, who took attendees through the history of Bioconductor, from the earliest contribution on the Bioconductor support site to the project’s growth into the global community it is today. The presentation highlighted how the project has evolved over the past 25 years and its impact on computational biology. The celebrations continued at the conference dinner, where attendees marked the occasion with a special anniversary cake. During the evening, Levi Waldron, one of Bioconductor’s Principal Investigators, shared a personal reflection on his journey with Bioconductor, from first encountering the project through his collaborations with Martin Morgan to becoming part of its leadership.

Levi Waldron shares a personal reflection on his journey with Bioconductor during the 25th anniversary celebrations.

Infrastructure and tools

Zulip

EuroBioC2026 continued to use Zulip as its primary communication platform. A dedicated conference channel, along with a separate hackathon channel, organised into topic-based threads, served as a central location for announcements, technical support, social interactions, and discussions before, during, and after the event. The threaded conversation model made it easier to follow discussions and kept participants connected throughout the conference.

Sticker Hexwall

The sticker hexwall returned for EuroBioC2026 following its successful introduction in 2025. The display showcased Bioconductor package stickers contributed by Bioconductor community members and served as a visual representation of the diversity of software projects within the ecosystem.

The hexwall quickly became a popular gathering point and photo location throughout the conference.

The hexwall at EuroBioC2026.

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

Le 16 juin 2026, Arthur Mensch, le patron emblématique de Mistral AI, a publié sur LinkedIn un message d’une rare solennité. Le discours officiel y célèbre, comme toujours, l’indépendance technologique européenne, la « souveraineté » et la mise à disposition imminente de nouveaux poids ouverts pour l’été.

Pourtant, sous le vernis des relations publiques et de l’alignement politique national, se glisse un aveu de taille :

« Aujourd’hui, nous ne possédons pas encore les meilleurs modèles de langage, mais nous avons constamment réduit cet écart. »

Cette phrase est une clé de lecture fondamentale. Elle signe la fin de l’illusion. L’écart technologique sur le raisonnement général pur face aux géants américains ne se comble pas par la simple élégance de l’algorithme ; il se creuse sous le poids des gigawatts de serveurs et des milliards de dollars de capitalisation. Pour survivre face à la Silicon Valley, le champion européen de l’intelligence artificielle vient d’opérer un pivot stratégique discret, mais d’une envergure historique : la Palantirisation.

En choisissant de dépêcher des troupes d’ingénieurs d’élite directement chez ses clients pour faire plier ses modèles à leurs processus métiers complexes, Mistral AI abandonne discrètement le rêve de l’éditeur logiciel pur à marges infinies pour embrasser la réalité humaine, physique et réglementaire de l’économie européenne.

1. La capitulation du « Pure Software »

Pour mesurer l’ampleur de cette mue, il faut se souvenir de la promesse de départ. En juin 2023, la note stratégique d’amorçage de Mistral AI — un document fondateur de sept pages rédigé par des transfuges de Meta et de Google DeepMind — esquissait un modèle d’affaires à faire saliver les plus grands fonds de capital-risque de la planète. La startup promettait une stratégie « Open Core » articulée autour d’un effet de réseau gratuit (des modèles ouverts téléchargeables par tous) redirigeant vers leur plateforme payante d’API logicielles.

C’était la promesse du modèle SaaS classique appliqué à l’ère cognitive : un coût marginal d’inférence décroissant, des équipes ultra-légères, pas d’infrastructure physique à gérer, et des marges brutes de près de 80 %. C’est sur cette thèse de scalabilité pure que les investisseurs ont valorisé Mistral à hauteur de 12 milliards de dollars en moins de trois ans.

Mais l’IA en entreprise a une propriété physique tenace : elle ne s’auto-installe pas.

Brancher un modèle de langage générique sur un système informatique d’entreprise hérité des années 1990 (legacy), s’assurer que les données ne fuitent pas chez un sous-traitant cloud, nettoyer les bases de données sémantiques et orchestrer des flottes d’agents fiables demande un effort d’ingénierie humaine titanesque. L’IA n’est pas un logiciel fluide ; c’est une infrastructure lourde.

La réponse d’Arthur Mensch ? L’introduction officielle du modèle FDE :

« Nous les aidons avec des pro-services (FDE est le nom chic), car c’est essentiel pour assurer la réussite de nos clients. »

L’acronyme FDE — Forward Deployed Engineers — est le secret de fabrication jalousement gardé d’un autre géant de la tech américaine : Palantir Technologies. Ce modèle consiste à détacher des ingénieurs d’élite directement dans les bureaux du client (que ce soit une banque d’affaires, un constructeur aéronautique ou un ministère de la Défense) pour câbler le logiciel aux sources de données réelles.

Pour Palantir, ce modèle d’affaires a longtemps été boudé par Wall Street, qui y voyait une vulgaire activité de « société de services » à forte intensité de main-d’œuvre et à faible scalabilité. Mais aujourd’hui, c’est le seul qui fonctionne en Europe. En adoptant les FDE, Mistral admet que pour vendre de l’IA sur le Vieux Continent, il ne suffit pas de mettre à disposition une clé d’API. Il faut envoyer des hommes pour construire l’ouvrage d’art.

2. Le poids de la matière et la revente de kilowattheures

Le second axe de cette transformation est d’ordre physique et financier. Au cours du premier trimestre 2026, Mistral AI a contracté un emprunt d’un montant inédit de 830 millions de dollars auprès de la BNP Paribas. Le collatéral de cette dette ? Pas des actions, mais de la matière brute : 13 800 puces de calcul de dernière génération (les GPU Nvidia GB300).

Pour une startup logicielle, se charger d’une telle dette d’infrastructure est un pari extrêmement risqué. Si les puces dorment ou si le prix d’inférence mondial s’effondre sous le coup du dumping des tokens des géants de l’open-source (comme le chinois DeepSeek), l’entreprise risque l’asphyxie financière.

Pour honorer ses créances et générer du cash-flow immédiat, Mistral a dû descendre dans la cave. Elle est passée de créatrice de modèles à fournisseur d’électricité cognitive :

« Suite à notre investissement dans l’infrastructure de calcul, nous proposons également des services cloud d’IA hébergés. »

L’entreprise utilise désormais ses serveurs excédentaires pour faire de la revente de puissance de calcul pure (GPU-as-a-Service), louant ses puces à de grands groupes industriels européens comme ASML, Ericsson ou l’Agence Spatiale Européenne pour qu’ils y pré-entraînent leurs propres réseaux via sa plateforme de calcul Forge.

C’est une intégration verticale inversée remarquable. Pour financer ses modèles de pointe, Mistral doit se comporter comme un gestionnaire de réseau électrique, vendant du kilowattheure de calcul en même temps que du token d’intelligence.

3. La souveraineté d’orchestration : l’interdépendance du silicium

Le discours d’Arthur Mensch insiste lourdement sur un niveau de service et de sécurité « totalement découplé des fournisseurs américains ». C’est le cœur du positionnement marketing de la startup : offrir aux acteurs de l’industrie, de la santé et de l’État un environnement de confiance face à l’hégémonie de Microsoft et de Google.

Pour certains observateurs, ce grand écart entre un discours d’indépendance et la réalité matérielle (des puces Nvidia GB300 conçues en Californie et gravées à Taïwan par TSMC, une structure de capital à forte composante anglo-saxonne) relève de la contradiction.

C’est oublier une vérité fondamentale que la physique et la géopolitique des semi-conducteurs nous imposent : dans l’industrie de pointe, l’autarcie à 100 % est une chimère.

La fabrication d’une seule puce de calcul IA est le produit d’une interdépendance radicale et d’une symbiose de monopoles nationaux. Les États-Unis conçoivent le design logiciel, le Japon détient la chimie de pointe et le silicium ultra-pur, l’Allemagne produit les optiques de Zeiss et les lasers de Trumpf, et les Pays-Bas, via le monopole absolu d’ASML, assemblent les machines de lithographie EUV. Et même ces machines néerlandaises dépendent de brevets et de technologies d’origine américaine (comme les sources de lumière Cymer).

Dans ce réseau à clés d’accès dispersées, exiger de l’Europe ou de Mistral qu’ils possèdent la chaîne de bout en bout est un contresens. La souveraineté réelle ne réside pas dans l’autarcie, mais dans le contrôle de verrous stratégiques spécifiques.

C’est ce qu’il convient d’appeler la souveraineté d’orchestration.

En se positionnant sur l’hébergement physique sur site (on-premise via Vibe) et l’intégration métier sur mesure par des ingénieurs FDE, Mistral ne prétend pas fabriquer du silicium souverain. Elle sécurise la couche finale et la plus critique de la pile : celle qui touche aux données d’affaires, au savoir tacite et aux secrets industriels de ses clients. Vous louez peut-être des puces soumises aux législations d’exportation américaines, mais vous gardez le contrôle total sur l’intelligence décisionnelle et la compliance de votre organisation. C’est un levier de confiance inestimable, et le seul pragmatique pour les entreprises régulées et les États européens (comme le prouve l’accord de défense signé en janvier 2026 avec le Ministère des Armées).

Conclusion : Un pivot pragmatique et risqué

Le pivot « Palantir » de Mistral AI n’est pas un aveu de défaite. C’est un éclair de génie tactique et de pragmatisme européen.

En acceptant la réalité des contraintes physiques de la course au calcul brut généraliste, Mistral choisit de s’ancrer là où la valeur réelle et durable se capture sur notre continent : dans la plomberie des processus industriels critiques, l’orchestration locale d’agents et la garantie d’une conformité de terrain. L’IA d’action ne se gagnera pas uniquement dans des laboratoires de recherche simulant une intelligence pure à coups de téraflops, mais sur le terrain, en codant auprès des systèmes legacy de l’économie réelle.

Cependant, ce choix opérationnel a un coût. L’ingénierie de solutions et les ingénieurs déployés (FDE) pèsent lourdement sur la marge brute de l’entreprise et sont complexes à faire monter à l’échelle par rapport à un pur modèle de distribution d’API logicielle. Si la baisse globale des prix de l’inférence se poursuit, Mistral devra rapidement stabiliser ses revenus récurrents issus de sa suite logicielle d’orchestration pour amortir le remboursement de sa dette d’infrastructure de 830 millions de dollars.

L’Europe n’aura peut-être pas son OpenAI scalable aux marges infinies d’un monopole d’accès d’API. Mais si le pari d’Arthur Mensch réussit, elle aura son Palantir de l’IA : un acteur capable de mailler la puissance de calcul physique à la haute précision de l’intégration industrielle. C’est le triomphe de la lucidité d’usage sur le fantasme de l’autarcie.

Sources et Références

• [1] Strategical Seed Memo of Mistral AI (juin 2023) : “Mistral AI: generative AI at European scale”, décrivant le business model de plateforme d’API logicielle Open Core à marges brutes scalables de type SaaS.

• [2] Benchmarks de Raisonnement d’Entreprise 2025-2026 : Évaluations d’inférence documentant l’écart persistant de 12 à 18 mois entre les modèles open-weights de Mistral AI et les capacités de raisonnement logique de GPT-4o, o1 ou Claude 3.5 Sonnet.

• [3] Palantir Technologies S-1 Registration Statement (2020) : Document d’introduction en bourse décrivant le modèle opérationnel des Forward Deployed Engineers (FDE).

• [4] Financement d’infrastructure BNP Paribas – Mistral AI (Q1 2026) : Structuration du crédit à hauteur de 830 millions de dollars garanti sur collatéral d’actifs physiques (13 800 GPU Nvidia GB300).

• [5] Accord-cadre Ministère des Armées – Mistral AI (janvier 2026) : Protocole de déploiement souverain de modèles de langage sur site pour les forces armées françaises.

A Likert scale (pronounced LICK-ert, not “LIKE-ert”) is a psychometric rating scale used in surveys and questionnaires to measure attitudes, opinions, and perceptions. Named after American social psychologist Rensis Likert, who developed it in 1932, it remains the most widely used approach to scaling responses in survey research today.

Likert Scale vs. Likert Item: An Important Distinction

These two terms are often used interchangeably — that is technically incorrect and matters for your analysis.

• A Likert item is a single statement with a rated response (e.g., “I am satisfied with the service: Strongly Disagree → Strongly Agree”).
• A Likert scale is the sum or average of several related Likert items designed to measure a single construct.

Treating a single item as a complete “scale” is one of the most common errors in survey design. If you have only one question, you have a Likert item, not a Likert scale — and the appropriate statistical treatment differs.

Types of Likert Scale Response Options

Likert scales are not limited to measuring agreement. Depending on your research objective, you can measure frequency, importance, quality, or likelihood using the same format. The table below shows the most common response option sets:

Likert Scale Formats: 4, 5, 6, and 7 Points Compared

Choosing the right number of response points affects the precision of your data and the cognitive load on respondents. Here is how each format differs in practice:

5-Point Likert Scale (Most Common)

The 5-point scale is the default choice in most survey research because it balances nuance with simplicity. Example:

“The quality of food at XYZ Restaurant is excellent.”

1. Strongly Disagree
2. Disagree
3. Neither Agree nor Disagree
4. Agree
5. Strongly Agree

4-Point Likert Scale (Forced Choice)

Removing the neutral midpoint forces respondents to take a position. Use this when fence-sitting would undermine your research objective — for example, when measuring purchase intent or policy support where “no opinion” is not useful data.

1. Strongly Disagree
2. Disagree
3. Agree
4. Strongly Agree

6-Point Likert Scale

Like the 4-point, this eliminates the neutral option while providing more granularity. Useful in employee satisfaction or consumer preference research where a clearer lean is needed.

1. Strongly Disagree
2. Disagree
3. Slightly Disagree
4. Slightly Agree
5. Agree
6. Strongly Agree

7-Point Likert Scale

The 7-point scale is preferred in academic and psychological research where capturing subtle differences in attitude matters. It improves statistical reliability but requires more careful labeling.

1. Strongly Disagree
2. Moderately Disagree
3. Slightly Disagree
4. Neither Agree nor Disagree
5. Slightly Agree
6. Moderately Agree
7. Strongly Agree

When to Use a Likert Scale

A Likert scale is the right tool when you need to measure characteristics that have no objective measurement — attitudes, opinions, satisfaction levels, or perceived likelihood. It is not appropriate when:

• A simple yes/no question would fully answer your research question
• You are measuring factual behaviors (e.g., “How many times per week do you exercise?” — use a numerical input instead)
• Respondents lack sufficient knowledge of the topic to have a genuine opinion

Use a Likert scale when you need to distinguish between degrees of agreement, not just direction. The difference between “Agree” and “Strongly Agree” often carries meaningful information in customer satisfaction and employee engagement research.

How to Design an Effective Likert Scale

Write Clear, Single-Focus Statements

Each item must address exactly one idea. A statement like “The service was fast and the staff were friendly” is a double-barreled item — the respondent may agree with one half and disagree with the other, making their response uninterpretable.

Balance Your Scale

A well-designed Likert scale includes an equal number of positively and negatively worded items. This counteracts acquiescence bias — the tendency of some respondents to agree with statements regardless of content. If all your items are positive, respondents who habitually agree will appear more satisfied than they actually are.

Avoid Leading Language

Avoid adverbs like “very,” “extremely,” or “always” inside the item statement itself. “This website is extremely fast” will yield fewer “Strongly Agree” responses than “This website is fast,” not because respondents think differently but because the bar is higher.

Keep the Scale Consistent Throughout the Survey

Switching between a 5-point and 7-point scale in the same questionnaire forces respondents to mentally reset and increases error rates. Choose one format and use it throughout.

Likert Scale Response Bias: What Can Distort Your Data

Understanding bias is not optional for anyone analyzing Likert data — it directly affects whether your conclusions are valid.

How to Analyze Likert Scale Data

Single Item vs. Multi-Item Scale: Different Rules Apply

This is the most commonly misunderstood part of Likert analysis. A single Likert item produces ordinal data — the intervals between response options are not guaranteed to be equal. Calculating a mean on ordinal data is statistically questionable. For a single item, use:

• Median as your measure of central tendency
• Frequency tables and percentages for distribution
• Chi-square tests or Mann-Whitney U for group comparisons

A full Likert scale (summed or averaged across multiple items) behaves more like interval data, especially with 5+ items and a reasonable sample size. In this case, parametric statistics become more defensible:

• Mean and standard deviation for descriptive summaries
• Cronbach’s alpha (α) to test internal consistency — aim for α > 0.7
• t-tests or ANOVA for group comparisons
• Spearman correlation for relationships between Likert scores and other variables

Cronbach’s Alpha: Checking if Your Scale Holds Together

If you are using multiple Likert items to measure the same construct, run Cronbach’s alpha before reporting results. An alpha above 0.8 indicates strong internal consistency. Values between 0.7 and 0.8 are acceptable. Below 0.7 suggests your items are not measuring the same thing — revise or remove items with low item-total correlations.

Likert Scale Examples Across Research Domains

Customer satisfaction survey:

“How satisfied are you with the cleanliness of our facilities?”

• Very Dissatisfied
• Dissatisfied
• Neither Satisfied nor Dissatisfied
• Satisfied
• Very Satisfied

Employee engagement survey:

“To what extent do you agree: ‘The new company policy enhances employee productivity’?”

• Strongly Disagree
• Disagree
• Neither Agree nor Disagree
• Agree
• Strongly Agree

Online UX research:

“Rate your agreement: ‘The online shopping experience was user-friendly and intuitive.'”

• Strongly Disagree
• Disagree
• Neither Agree nor Disagree
• Agree
• Strongly Agree

Advantages and Disadvantages of Likert Scales

Conclusion

A Likert scale is one of the most versatile and reliable tools in survey research — when used correctly. The key decisions are choosing the right number of response points for your research goal, writing items that are balanced and unambiguous, and applying the correct statistical method depending on whether you are working with a single item or a multi-item scale. Whether you are measuring customer satisfaction, employee engagement, student attitudes, or any other opinion-based construct, the principles remain the same: clarity in item wording, consistency in format, and honesty about what ordinal data can and cannot tell you.

Frequently Asked Questions

Q: What is the difference between a Likert scale and a Likert item?

A: A Likert item is a single rated statement. A Likert scale is the aggregate of multiple related items. The distinction matters for analysis: single items should use median and nonparametric tests; full scales can use mean and parametric tests.

Q: How do you pronounce Likert?

A: The correct pronunciation is “LICK-ert,” not “LIKE-ert.” It is named after Rensis Likert, who created the scale in 1932.

Q: Should I use a 5-point or 7-point Likert scale?

A: For general surveys and applied research, a 5-point scale is easier for respondents and produces reliable results. For academic or psychological research where detecting subtle attitude differences matters, a 7-point scale offers better statistical discrimination. Research comparing 5-point and 7-point scales finds that both produce similar mean scores once rescaled — so the choice depends more on respondent context than statistical superiority.

Q: Can you calculate the mean from Likert scale data?

A: For a single Likert item, technically no — the data is ordinal, so the median is more appropriate. For a complete Likert scale (multiple items summed), calculating the mean is widely practiced and generally acceptable, especially with a sample size above 30 and if the data distribution is approximately normal.

Q: What is acquiescence bias in Likert scales?

A: Acquiescence bias is the tendency of some respondents to agree with statements regardless of content. It is reduced by including both positively and negatively worded items in your scale, so that habitual agreement on one item is balanced by habitual agreement on an item that pulls in the opposite direction.

Q: Are Likert scale questions suitable for all types of research?

A: Likert scales work well in social sciences, market research, psychology, education research, healthcare, and UX research. They are not appropriate when you need objective behavioral counts or factual data — use open-ended questions or numerical inputs for those cases.

Q: Is it necessary to include a neutral response option in a Likert scale?

A: No. Including a neutral option (odd-point scale) allows genuinely ambivalent respondents to express that accurately. Removing it (even-point scale) forces a directional choice, which can reduce central tendency bias but may frustrate respondents who truly have no strong view. Choose based on whether neutrality is meaningful in your research context.

The Rousseeuw Prize honors five pioneering developers for nearly three decades of unpaid work building R, the foundational open-source computing language behind artificial intelligence, healthcare, and economic decision-making.

• The $1 million Rousseeuw Prize for Statistics recognizes three decades of foundational work that transformed how statistical methods are developed, validated, and shared globally.
• R, the open-source statistical computing language, underpins modern AI development, pharmaceutical research, financial modeling, and global scientific analysis.
• Used by organizations including the U.S. Food and Drug Administration, major pharmaceutical companies, and global central banks, R has become the trusted infrastructure for high-stakes analysis because it is stable, auditable, and reproducible.

NEW YORK – June 17, 2026 — Five members of the R Core Team have been awarded the prestigious Rousseeuw Prize for Statistics for their decades of work building and maintaining the R Project, “R,” a free and open-source statistical computing language used across global research institutions, healthcare systems, financial organizations, and technology companies. The Rousseeuw Prize is an international award recognizing major contributions to statistical research. 

The 2026 Rousseeuw Prize honorees are:

• Brian D. Ripley, emeritus professor at the University of Oxford
• Martin Maechler, emeritus professor at ETH Zurich
• Kurt Hornik, department chair at WU Vienna University of Economics and Business
• Peter Dalgaard, professor at Copenhagen Business School
• Luke Tierney, professor at the University of Iowa

The five laureates receive half of the prize money because they are deemed to have made the longest sustained contributions to the R project; the other half of the prize is shared among the many others who have been active on the R Core Team.

Together, the R Project volunteers have spent the last 27 years and a collective 28,000 coding hours on R, developing an open-source programming language and software environment that transformed statistics from a proprietary corporate tool into a global public good. The software is relied upon by organizations including the U.S. Food and Drug Administration, pharmaceutical companies, and central banks such as the European Central Bank and the Bank of England.

The award recognizes the team’s role in making advanced statistical tools widely accessible. By keeping R free and open-source under the GNU General Public License, the R Core Team removed many of the financial barriers that have historically limited access to advanced analytics software. Due to this increased accessibility, hundreds of thousands of users including researchers, students, hospitals, public health organizations, and governments around the world are able to utilize the same statistical tools regardless of institutional resources. In addition, they use R to share transcripts of their data analyses, allowing one user’s workflows to power other users data analyses everywhere around the world. The frictionless spread of these transcripts has powered countless educational data science projects globally and hundreds of course textbooks at the PhD and Master’s level. In a recent twist, it’s not only humans who use R: AI data analyst `agents’ have been learning from the massive volume of published R transcripts and are now able to assist with many everyday data analysis tasks. 

“Long before AI became a global conversation, the R Core Team was building the statistical infrastructure that made today’s data-driven world possible,” said Stanford University statistics professor and leading statistician David Donoho, PhD. “This team’s stewardship of R created an open and trusted foundation for research across disciplines and continents. Few innovations have had such a profound effect on how knowledge is produced, shared, and validated in the modern era.”

Named after Professor Peter Rousseeuw, a pioneering Belgian statistician known for his foundational work in robust statistics and data analysis, the Rousseeuw Prize for Statistics recognizes innovations that have transformed the understanding and application of data for the benefit of society. Past laureates include internationally renowned statisticians and researchers whose work has advanced fields ranging from epidemiology and artificial intelligence to public policy and scientific discovery.

For more information, visit https://www.rousseeuwprize.org/.

###

Media Contact:

rousseeuwprize@ampublicrelations.com

{talib} is a new R package built on TA-Lib, which is now available on CRAN. The R-package is targeted at individuals and, perhaps, institutions who, in some form or the other, interacts with the financial markets using technical analysis.

The library is built with minimal dependencies for long-term stability and freedom in mind. All functions are built around data.frame– and matrix-classes which are portable to all other data-containers with minimal effort.

Everything in the library is built ‘bottom-up’ for maximum speed and memory efficiency. Each indicator interacts directly with R’s C API via .Call().

In this blog post I will give a brief introduction to the charting interface which is built to mimick the behaviour of base R’s plotting API.

talib::chart(
  talib::BTC
)

str(formals(talib::chart))
#> Dotted pair list of 5
#>  $ x    : symbol 
#>  $ type : chr "candlestick"
#>  $ idx  : NULL
#>  $ title: symbol 
#>  $ ...  : symbol

talib::set_theme("hawks_and_doves")

talib::chart(
  talib::BTC
)

{
  talib::chart(talib::BTC)
  talib::indicator(talib::SMA, n = 7)
  talib::indicator(talib::SMA, n = 14)
  talib::indicator(talib::SMA, n = 21)
  talib::indicator(talib::SMA, n = 28)
  talib::indicator(talib::MACD)
  talib::indicator(talib::trading_volume)
}

install.packages("talib")

install.packages(
  "talib",
  type = "source",
  configure.args = "-O3 -march=native"
)

shiny.webawesome brings Web Awesome to R Shiny through generated wrappers, reactive bindings, and a bundled runtime. It aims for complete component support while staying close enough to upstream that the Web Awesome docs and examples are directly useful in everyday package use.

CRAN | R-universe | Package website | Source repository

Background

shiny.webawesome started from a perceived gap: Shiny would benefit from a UI library that feels modern, visually polished, and broad enough to support a full app coherently. Web Awesome was a strong fit because it combines rich components, layout and styling utilities, and detailed upstream documentation with a standards-based web-components structure that is straightforward to track from R. That makes it easier for the package to stay close to upstream while still fitting naturally into Shiny.

The Whole Game

Here’s a screenshot of a simple, complete example app using shiny.webawesome. The full live app and code are available in an article at: https://mbanand.github.io/ghpages/announcement/..

This example showcases many of the facilities available in the package:

• a visually rich component library
• direct use of Web Awesome layout utilities such as wa-stack, wa-cluster, wa-gap-*, and wa-align-* classes
• styling through Web Awesome design tokens and classes such as --wa-color-*, --wa-font-*, and wa-body-*
• reactive Shiny input bindings
• helpers for calling methods on HTML elements, setting properties, and injecting simple JavaScript snippets

Design Philosophy

shiny.webawesome is designed to stay close to upstream Web Awesome. Most component wrappers are generated from Web Awesome metadata, which helps preserve upstream names, structure, and behavior while translating the interface into normal R conventions such as snake_case.

That close alignment has a practical benefit: when you want deeper details, examples, or component-specific guidance, you can usually go straight to the upstream Web Awesome documentation and apply what you find directly in shiny.webawesome. The package currently supports all Web Awesome components, so the upstream docs are a practical reference for day-to-day use.

To support the server-client model of Shiny, the package adds a small set of page and layout helpers, curated reactive bindings, and a narrow command layer for cases where browser-side interaction goes beyond the generated wrappers.

The result is a package with a clear default path. Use generated wrappers for ordinary UI, use bindings for meaningful reactive state, and reach for commands or small JavaScript glue when the app needs them.

Shiny Bindings

shiny.webawesome does not forward every browser event and every detail of component telemetry into Shiny. Much component state and interaction detail is better handled locally in the browser rather than turned into server messages. Consequently, the package exposes only a curated set of Shiny bindings that fit Shiny’s reactive model, with an emphasis on meaningful committed state rather than low-level browser event streams.

In the most common case, a binding publishes a durable semantic value. A select reports its current value, a dialog can report whether it is open, and a tree can report the currently selected item ids. The key idea is that Shiny receives the state the app actually cares about, not the raw event name that happened to produce it.

Some components are better treated as actions than values. A button is the clearest example: in Shiny, it behaves like a Shiny action input, with each click producing a new input event. A small number of components need both action semantics and a separate value. A dropdown, for example, may need to trigger reactivity on every choice, including repeated selections of the same item, while also exposing the latest selected value.

This design keeps reactive messaging to the server smaller, clearer, and easier to reason about. If an interaction belongs naturally in Shiny’s input model, shiny.webawesome will expose it as a binding. If it is more naturally a browser-side concern, it is usually a better fit for the command layer or a small amount of JavaScript glue.

For the full binding categories, semantics, and examples, see the package article: Shiny Bindings.

Command API

shiny.webawesome covers the most common interaction patterns through generated wrappers, Shiny bindings, and update helpers. But sometimes an app still needs to reach into a live browser element directly: set a property, call a method, or add a small browser-local JavaScript snippet.

For those cases, the package provides a narrow command API. The two main server-side helpers are wa_set_property() and wa_call_method(). They let Shiny code send one-way commands to a browser element identified by id, either by assigning a value to a live property or invoking a browser-side method.

If a component already has a binding or update helper, that should usually remain the first choice. The command layer is for the cases that fall just outside those built-in paths, where the simplest solution is still to tell the existing browser component to do one specific thing.

The package also includes wa_js() for a different kind of job: small, app-local JavaScript glue. That is useful when the missing piece is browser-side logic such as listening for an event, reading live component state, or publishing a derived value back to Shiny with Shiny.setInputValue().

For more detail and examples, see the package article: Command API.

Conclusion

shiny.webawesome brings a visually rich component library into Shiny while staying close to upstream Web Awesome. That combination gives polished components, useful layout and styling utilities, and a workflow where upstream documentation and examples remain directly relevant throughout app development.

For more examples, longer articles, and full reference material, see the package website: shiny-webawesome.org.

Introduction

If you’re an R user, you’ve probably experienced these moments:

• You’re writing code and forgot the exact syntax for a function
• Your code throws an error and you’re staring at a confusing error message
• You have a block of code but want to understand what it does in plain English
• You want to chat with an AI assistant about your data analysis, but don’t want to leave RStudio

llmcoder is an RStudio addin that solves all of these problems by integrating Large Language Model (LLM) assistance directly into your RStudio workflow, and more importantly, it’s FREE!

In this post, I’ll show you how llmcoder can speed up your R coding and make your workflow smoother.

Watch a quick demo of llmcoder in action:

https://youtu.be/SRzjaURbKCw

Installation

You can install llmcoder from GitHub:

# Install remotes if you haven't already
install.packages("remotes")

# Install llmcoder
remotes::install_github("ShiyangZheng/llmcoder")

Load the package:

library(llmcoder)

Feature 1: Generate R Code from Inline Comments

Ever wish you could just type what you want in plain English and get R code instantly?

How to use:

1. Type a comment describing what you want
2. Place your cursor on that line
3. Use the Addins menu and select “Generate Code from Comment”

Example:

# Load the mtcars dataset and create a scatter plot of mpg vs wt, colored by number of cylinders

After running the addin, the comment is replaced with:

library(ggplot2)
data(mtcars)
ggplot(mtcars, aes(x = wt, y = mpg, color = factor(cyl))) +
  geom_point(size = 3, alpha = 0.8) +
  labs(
    title = "Fuel Efficiency vs Weight by Cylinder Count",
    x = "Weight (1000 lbs)",
    y = "Miles per Gallon",
    color = "Cylinders"
  ) +
  theme_minimal()

No more switching to ChatGPT or copying code from Stack Overflow!

Feature 2: Fix Console Errors with LLM Assistance

We’ve all been there – a cryptic error message and you’re not sure what went wrong.

How to use:

1. Run code that produces an error
2. The error appears in the console
3. Use the Addins menu and select “Fix Error with LLM”

Example:

library(dplyr)
data %>%
  filter(cyl == 4) %>%
  summary()
# Error: object 'data' not found

llmcoder captures the error and sends it to the LLM, which returns an explanation and suggests:

mtcars %>% filter(cyl == 4) %>% summary()

Feature 3: Explain Selected Code in Plain English

Sometimes you inherit code from a colleague or find a Stack Overflow answer and want to understand what it does.

How to use:

1. Select a block of code in the editor
2. Use the Addins menu and select “Explain Code”

Example:

mtcars %>%
  group_by(cyl) %>%
  summarize(
    mean_mpg = mean(mpg, na.rm = TRUE),
    sd_mpg = sd(mpg, na.rm = TRUE),
    count = n()
  ) %>%
  arrange(desc(mean_mpg))

llmcoder returns:

1. Takes the built-in mtcars dataset
2. Groups the data by the number of cylinders (cyl)
3. Calculates the mean and standard deviation of miles per gallon (mpg) for each group
4. Arranges the results in descending order of mean fuel efficiency

Feature 4: Multi-Turn Chat Panel with Session Context

This is the flagship feature. llmcoder includes a Chat Panel that understands your current R session.

How to open: Use the Addins menu and select “Open Chat Panel”

What makes it special?

The Chat Panel is session-aware:

• It knows which packages you have loaded
• It knows what objects are in your global environment
• It can read the contents of your current script
• It has access to your recent console history

Example conversation:

You: What’s the correlation between mpg and wt in mtcars?

AI: The correlation between mpg and wt in the mtcars dataset is -0.87, indicating a strong negative relationship. As weight increases, fuel efficiency decreases.

cor(mtcars$mpg, mtcars$wt, use = "complete.obs")

Want to see the Chat Panel in action? Watch this demo:
https://youtu.be/zP-RuCN3q14

Supported LLM Providers

llmcoder supports multiple LLM providers – you can choose the one that works best for you:

Privacy note: If you use Ollama, all processing happens locally on your machine. No data is sent to external servers.

Customization: Choose Your Prompt Style

The Chat Panel allows you to select different prompt styles:

• General Assistant: Best for general questions
• R Code Helper: Focuses on writing clean, idiomatic R code
• Statistics Advisor: Helps with statistical concepts and test selection
• Research (Psycho): Tailored for psycholinguistics researchers

Why llmcoder?

There are many AI coding assistants out there (Copilot, Cursor, etc.), so why llmcoder?

1. Native RStudio integration: No need to switch to another app or browser tab
2. Session-aware: The LLM knows what you’re working on
3. Multiple LLM providers: Choose the one you prefer (or use a local model for privacy)
4. Open source: MIT license, free to use and modify
5. Designed for R users: Not a generic coding assistant – it understands R-specific workflows

Call to Action

Ready to try llmcoder?

remotes::install_github("ShiyangZheng/llmcoder")

GitHub: https://github.com/ShiyangZheng/llmcoder

If you encounter any bugs or have feature requests, please file an issue: https://github.com/ShiyangZheng/llmcoder/issues

Star the repo if you find it useful!

About the Author

Shiyang Zheng is a PhD student in Psycholinguistics at the University of Nottingham. His research focuses on idiom acquisition and computational modeling. He built llmcoder to make R coding easier for himself and the R community.

• GitHub: @ShiyangZheng
• Academic website: shiyangzheng.top
• ORCID: 0000-0003-0511-4683

Clinical data are rarely as clean, compact or convenient as the examples we often use when teaching statistics or R. Real hospital datasets are usually distributed across several tables, include missing values, contain repeated structures, and require careful documentation before they can be reused.

The new R package DIVINE is interesting precisely because it brings that reality into the R ecosystem in an accessible way.  Available on CRAN, DIVINE provides a curated collection of datasets from a multicentre cohort of hospitalized COVID-19 patients in the south metropolitan area of Barcelona. The package is accompanied by a recent publication in Scientific Data, which describes the database, its structure, data collection process and potential reuse for clinical epidemiology, teaching and methodological research.

A clinical dataset packaged for R

The package includes 14 datasets covering different clinical domains, such as demographics, comorbidities, symptoms, vital signs, severity scores, ICU information, treatments, complications, vaccination and end-of-follow-up data.

This relational structure is one of the most valuable aspects of the package. Instead of providing a single pre-merged analysis file, DIVINE preserves the logic of a real clinical database, where information is distributed across several linked tables. This makes it especially useful for applied teaching and for demonstrating realistic data-management workflows in R.

For example:

install.packages("DIVINE")
library(DIVINE)

data(package = "DIVINE")

The datasets can then be loaded in the usual way:

data("demographic")
data("vital_signs")
data("scores")

The common identifiers allow users to combine information across tables and build analysis datasets depending on the research question.

More than a data package

Although the datasets are the main contribution, DIVINE also includes helper functions for common epidemiological data workflows. These include:

data_overview()
multi_join()
stats_table()
multi_plot()
impute_missing()
export_data()

These functions are not intended to replace the broader R ecosystem, but they make the package easier to use in teaching, exploratory analysis and reproducible examples.

A minimal workflow might look like this:

library(DIVINE)

data("demographic")
data("vital_signs")
data("scores")

baseline <- multi_join(
  list(demographic, vital_signs, scores),
  key = c("record_id", "covid_wave", "center"),
  join_type = "left"
)

data_overview(baseline)

stats_table(
  baseline,
  vars = c("age", "sex"),
  by = "covid_wave",
  statistic_type = "median_iqr",
  pvalue = TRUE
)

This example already illustrates several important aspects of clinical data analysis: understanding table structure, joining related datasets, checking variables, and producing descriptive summaries.

Why it is useful for R users

For a specialised R audience, the value of DIVINE is not only that it provides COVID-19 data. Its main interest is that it offers a realistic, documented and reusable clinical database within a familiar R workflow.

The package may be useful for:

• teaching data management with relational clinical datasets;

• preparing examples for biostatistics or epidemiology courses;

• demonstrating descriptive clinical analyses;

• exploring missing data and variable availability;

• developing prognostic modelling examples;

• validating prediction models;

• creating reproducible workflows using real-world health data.

This makes DIVINE particularly attractive for applied biostatisticians, epidemiologists, clinical researchers and R instructors who want to move beyond toy datasets.

Snapshot testing is not about screenshots.

Most people meet it through UI regression tests: render a component, save a picture, fail the build when the picture changes. So the technique gets filed away as “the thing that compares images.” That is one use. But not the only one.

The mechanic underneath is general. Capture some output, save it to a file, and on every later run compare fresh output against the saved copy. The output can be a plot. It can also be console text, a log, a data frame, an error message, or a deeply nested list. Anything you can serialize, you can snapshot.

What makes it powerful is also what makes it dangerous: you are the test oracle. There is no expect_equal(result, 42) stating the answer up front. You accept the first snapshot because you read it and judged it correct. Get that review wrong, or skip it, and you have pinned a bug in place and called it a passing test.

In this post I want to walk through using snapshot testing for what it is good for, and the practices that make it efficient.

What snapshot testing actually is

In testthat’s third edition the entry points are expect_snapshot() and expect_snapshot_file(). The first run records output into a _snaps/ directory next to your tests, as a .md file named after the test file. Every run after that compares against what’s recorded. A mismatch fails the test and shows you a diff.

test_that("summary prints a one-line overview", {
  expect_snapshot(print(summary(1:10)))
})

The first time, testthat writes the printed output to _snaps/summary.md and the test passes (with a note that a new snapshot was recorded). From then on, that file is the expected value.

You reach for this in those situations:

• The output is large or tedious to assert field by field, but you can recognize whether it’s correct by looking at it.
• The output is impossible to express in code: a rendered plot, a rendered table, any image.
• The output is impractical to express in code a formatted CLI report, a full console transcript with its alignment. You can’t write expect_equal() for “the table is laid out correctly.” You can look at it and know.

Five practices that keep snapshots trustworthy

Snapshot suites rot in predictable ways: noise the diff engine flags as failures, snapshots nobody can review, tests that flake on a different machine, and (as with any other test) titles that say nothing.

Those practices prevent it.

1. Scope to exactly what proves the behavior

Capture the plot, not the page the plot lives on.

If you’re testing that a chart colors points correctly, snapshot the chart. Not the dashboard it’s embedded in, with its header, its sidebar, the current date in the corner, and a “last refreshed” timestamp. Every one of those is unrelated to the behavior under test, and every one is a reason for the comparison to fail when nothing you care about changed.

A snapshot’s diff engine is literal. It flags any difference. So don’t hand it differences that don’t matter. Scope the capture down to the smallest thing that demonstrates the behavior, and the only way the test can fail is if that behavior breaks.

Don’t give the diff engine extra reasons to make false positives.

2. Make snapshots human-readable

You are going to review these files by eye. So they have to be readable by eye.

Store snapshots as text: markdown, CSV, SVG, JSON. Never as a binary blob. I heard on the R Weekly podcast that some teams keep snapshots as .rds files, and I’d push back on that hard. A binary snapshot can’t be read in an editor, can’t be reviewed in a pull request, and can’t be diffed when it changes. It defeats the entire premise. The whole technique rests on a human being able to look at the recorded output and decide it’s right. You want to also help your code reviewers to do that, don’t hide the “truth” in a binary file. Especially when accepting the first snapshot as “the truth”; make it easy for your collaborators to read and judge the snapshot!

Don’t introduce extra points of friction. Keep it simple.

3. Remove nondeterminism, or filter what’s left

A snapshot that changes on every run is useless. Timestamps, random IDs, elapsed-time measurements, unordered query results. Any of these will make the file churn and train you to accept changes blindly.

Fix it at the source first. Inject the things that vary so the test controls them: pass a fixed clock instead of calling Sys.time(), set a seed, supply IDs rather than generating them. This is dependency injection, the same move that makes any code testable.

When you can’t remove the variation, filter it. testthat’s expect_snapshot() takes a transform argument: a function that cleans each line of output before it’s compared. Strip the timestamps, drop the spinner characters, normalize the paths. For data, impose a deterministic order before you serialize.

Don’t let snapshots change when there is no reason for them to change.

4. Stabilize platform differences

A rendered snapshot depends on more than your code. Fonts render differently on macOS and Linux. A new release of a plotting or formatting dependency shifts the output by a pixel or a label. R itself changes between versions. None of that is a regression, but a literal diff engine can’t tell, so a snapshot recorded on your laptop fails the moment it runs anywhere else.

Two tools handle this, and they work together.

A. Variants keep incompatible environments from overwriting each other. Both expect_snapshot() and expect_snapshot_file() take a variant argument. testthat stores each variant in its own subdirectory, _snaps/{variant}/, so the macOS render and the Linux render sit side by side instead of clobbering one another. Key the variant on whatever actually moves the output: the operating system, the R version, a specific dependency’s version, or a combination. You decide what relevant axes of variation are, and you key the snapshots to them. Maybe you want to support rendering on different platforms, and you want to support different versions of a plotting library. Then the variant should include both the platform and the library version.

variant = paste(platform_variant(), packageVersion("echarts4r"), sep = "-")

B. Let one platform generate the truth, and let the whole team use it. Variants solve the storage problem. They don’t solve the contribution problem: your developers are on different operating systems, and you don’t want each regenerated snapshot to depend on whose machine produced it. If your team works on Windows, macOS and Linux, you may not want to check-into the repository 3 slightly different copies of the same thing. Nominate a single canonical environment, your CI runner, and treat the snapshots it produces as authoritative.

When a snapshot test fails on GitHub Actions, the files that run produced are uploaded as build artifacts. testthat gives you a helper to pull them straight into your local checkout:

testthat::snapshot_download_gh(
  repository = "your-org/your-package",
  run_id = "47905180716"
)

This is a quite recent addition to testthat. Worth knowing.

You rarely have to look the call up. When snapshots fail inside an R CMD check job, testthat prints the exact snapshot_download_gh() line in the CI log, ready to copy. Run it, review the downloaded files the way you’d review any first snapshot, and commit them.

That turns snapshot testing into a team practice. A contributor on Windows can change a plot, open a pull request, and let CI render the canonical image. The reviewer accepts the snapshot CI produced, not one tied to a particular laptop. The truth comes from one place, and everyone contributes to it through the same door.

You’ll notice both options have their advantages and disadvantages. It’s up to you to decide which one fits your team and workflow better.

But now that you know your options you can test them out and see which one works best for you.

5. Name the test and the snapshot so they stand alone

A test title should state the precondition and the expected output.

The same holds for snapshot tests. Not “reporter works.” Something like “progress reporter shows survived mutants in summary.” The title is the first thing a reviewer reads when the snapshot changes.

But snapshot tests aren’t self contained.

The assertion of a snapshot test is a file. That means you read the test and then you need to open the file to understand what is really the expected outcome. But there is also another workflow: you might also browse the snapshots directory first and get a grasp of what the code is producing.

With expect_snapshot_file() you also name the snapshot file yourself. Use that. A file called scatterplot_colors_points_in_the_band.png tells you what it should contain before you even read the test itself. The filename and its content should tell the story on their own, without you having to dig up the test that produced them.

The rest of this post is five worked examples, each leaning on these five practices.

Example 1: plots, from simple to interactive

ggplot: snapshot the SVG, not a PNG

For ggplot, the right tool is vdiffr. It renders a plot to SVG, which is text, and snapshots that.

test_that("points below lower threshold are green, above upper are red, inside are yellow", {
  # Arrange
  data <- data.frame(
    x = seq_date(3),
    y = c(10, 20, 30)
  )

  # Act
  p <- threshold_plot(data, lower = 15, upper = 25)

  # Assert
  vdiffr::expect_doppelganger("threshold_plot_below_threshold_green_above_threshold_red_inside_threshold_yellow", p)
})

Two of our practices fall out of this for free. The snapshot is scoped to the plot object itself, not a Shiny page that embeds it. And it’s human-readable: the recorded .svg is text you can open, and vdiffr ships a Shiny app (vdiffr::manage_cases()) that shows the old and new render side by side when something changes. You review the picture, but the artifact under version control is inspectable text.

Here’s the actual plot that test captures:

And here are the first lines of the SVG vdiffr would record. This is the whole point of practice #2: the snapshot under version control is text you can read.

<?xml version="1.0" encoding="UTF-8"?>
<svg xmlns="http://www.w3.org/2000/svg" xmlns:xlink="http://www.w3.org/1999/xlink" width="450" viewBox="0 0 504 288">
<defs>
<g>
<g id="glyph-0-0">

Comparing text is easier than comparing pixels. If saving image to SVG is possible, you should always prefer it to a PNG. Not only you can view SVG both as text and the image, but also the diff engine can tell you exactly what changed in the markup, instead of just showing a pixel difference.

But there are plenty of cases when SVG isn’t an option.

htmlwidgets: when you’re forced back to pixels

vdiffr can’t render an htmlwidget or a Shiny tag list, because there’s no static SVG to produce. Here you fall back to rendering the thing to a real PNG and comparing images. That’s a harder problem, and it’s where scoping and determinism stop being nice-to-haves.

The shape of the helper: render to HTML, screenshot to PNG with webshot, then hand the PNG to expect_snapshot_file().

expect_plot <- function(x, name, ...) {
  UseMethod("expect_plot")
}

expect_plot.htmlwidget <- function(
  x,
  name,
  variant = shinytest2::platform_variant(),
  width = 992,
  height = 744
) {
  local_edition(3)
  html_temp <- fs::path(tempdir(), name, ext = "html")
  png_temp <- fs::path(tempdir(), name, ext = "png")
  on.exit(unlink(tempdir()))

  htmlwidgets::saveWidget(
    x,
    file = html_temp,
    selfcontained = FALSE
  )
  webshot::webshot(
    url = html_temp,
    file = png_temp,
    delay = 0.5,
    quiet = TRUE,
    vwidth = width,
    vheight = height
  )

  testthat::expect_snapshot_file(
    png_temp,
    name = fs::path(name, ext = "png"),
    variant = variant
  )
}

I’ve used this pattern across many projects and it always worked for me very well.

The threshold. A pixel-exact comparison of a rendered chart will fail on trivial, invisible differences in anti-aliasing. So the compare function allows a small per-pixel difference budget before it calls a mismatch. Locally (interactive()) the threshold is 0, because you want to see every change as you work. On CI it’s relaxed, because the CI renderer isn’t identical to your laptop.

The variant. Fonts render differently across operating systems, so a snapshot recorded on macOS will not match one produced on Linux. platform_variant() keys the snapshot to the platform, and expect_snapshot_file() keeps a separate recorded file per variant (_snaps/mac/..., _snaps/linux/...), so cross-platform rendering differences never masquerade as a regression. This is the storage half of practice #4; pair it with snapshot_download_gh() so CI generates the canonical files the whole team commits.

The test that uses it reads like any other, with a descriptive title and a named snapshot:

it("colors timepoints between thresholds", {
  # Arrange
  data <- data.frame(
    x = seq_date(3),
    y = c(10, 20, 30)
  )

  # Act
  plot <- threshold_plot(data, lower = 15, upper = 25)

  # Assert
  expect_plot(
    plot,
    name = "colors_timepoints_between_thresholds"
  )
})

Interactive plots: you can even snapshot an interaction

An htmlwidget is just HTML and JavaScript. That means you can drive it into a specific state and snapshot that. Here’s an echarts4r line chart where the behavior under test is the tooltip content.

A small helper dispatches the chart’s own “show tooltip” action when the widget loads, it simulates user interacting with the plot:

#' tests/testthat/setup-trigger_tooltip.R
trigger_tooltip <- function(x, series_index, data_index) {
  htmlwidgets::onRender(
    x,
    sprintf(
      "function(el, x, data) {
        const chart = echarts.getInstanceByDom(el);
        chart.dispatchAction({
          type: 'showTip',
          seriesIndex: %s,
          dataIndex: %s,
        });
      }",
      series_index,
      data_index
    )
  )
}

Then the test renders the chart, triggers the tooltip, and snapshots the result (notice that the # Act is triggering the tooltip, not creating the plot):

it("shows tooltip content from the specified tooltip column", {
  # Arrange
  data <- data.frame(
    date = as.Date(
      c("2020-01-01", "2020-02-01", "2020-03-01")
    ),
    value = c(1, 3, 2),
    tooltip = "TOOLTIP"
  )
  plot <- line_plot(
    data,
    date = "date",
    value = "value",
    tooltip = "tooltip"
  )

  # Act
  result <- plot %>%
    trigger_tooltip(series_index = 0, data_index = 1)

  # Assert
  expect_plot(result, "line_plot_tooltip")
})

The tooltip text is doing double duty. It’s the data the test checks, and it’s a human-readable marker that tells the reviewer exactly what to look for when accepting the snapshot. Here’s the PNG that snapshot captures:

Example 2: printed output (reporters, loggers, CLI)

Some objects exist to print.

Test reporters, loggers, CLI tools. Their whole job is to render formatted text to the console.

That text is the behavior, and expect_snapshot() captures it verbatim into a readable .md file.

I use this in muttest, a mutation testing package, to pin down exactly what the progress reporter prints. The test is small:

test_that("progress reporter shows all killed", {
  .with_example_dir("shipping/", {
    mutators <- list(operator(">", "<"))
    plan <- muttest_plan(mutators, fs::dir_ls("R"))
    .expect_snapshot(
      muttest(
        plan,
        reporter = ProgressMutationReporter$new(
          min_time = Inf,
          survived_detail = "none"
        )
      )
    )
  })
})

But a reporter’s output is full of nondeterminism: spinner frames, blank lines, and per-step timings like [0.3s] and Duration: 1.2s. Snapshot that raw and it fails on every run. The fix is practice #3: a transform that strips the noise before comparison. I wrap expect_snapshot() once, in setup.R, so every test in the suite gets the cleaning for free:

.expect_snapshot <- purrr::partial(
  testthat::expect_snapshot,
  transform = function(lines) {
    lines |>
      stringr::str_subset("^[\\|/\\-\\\\] \\|", negate = TRUE) |>
      stringr::str_subset("^$", negate = TRUE) |>
      stringr::str_remove_all("\\s\\[\\d+.\\d+s\\]") |>
      stringr::str_remove_all("Duration:\\s\\d+.\\d+\\ss") |>
      stringr::str_trim()
  }
)

The first two str_subset calls drop spinner lines and blanks. The two str_remove_all calls delete the timing fragments. What’s left is the stable, meaningful part of the output, and that’s what lands in the snapshot:

# progress reporter shows all killed

    Code
      muttest(plan, reporter = ProgressMutationReporter$new(min_time = Inf,
        survived_detail = "none"))
    Output
      i Mutation Testing
        |   K |   S |   E |   T |   % | Mutator  | File
      v |   1 |   0 |   0 |   1 | 100 | > → <    | shipping.R
      -- Results ---------------------------------------------------------------------
      [ KILLED 1 | SURVIVED 0 | ERRORS 0 | TOTAL 1 | SCORE 100.0% ]

This is a snapshot doing exactly what it should. The table alignment, the symbols, the score line. None of that is pleasant to assert by hand, but all of it is obviously correct (or obviously wrong) at a glance. The full reporter source, setup, and recorded snapshots are in the muttest repo: setup.R, the test, and the snapshot.

Example 3: data frames as CSV

Data frames are a classic snapshot candidate, and a classic way to get it wrong.

The tempting move is expect_snapshot(print(df)). Don’t. Printed data frames are truncated past a certain size, formatted to your console width, and shown in whatever order the rows happen to be in. You’re snapshotting the print method, not the data.

Write the data frame to CSV. CSV is text, it diffs cleanly, and it’s the obvious human-readable representation of tabular data.

I find snapshotting tables especially useful when you need signoff of business logic calculation from a business expert. Then instead of showing a table created in code you can print hand over the CSV or even a formatted markdown table for review. The expert can then sign off on the calculation without needing to read the code.

Following the same S3 pattern as expect_plot, here’s a custom expectation. The comparison is the part worth getting right, so it lives in its own named function:

compare_df <- function(old, new) {
  # Compare parsed data frames, not raw CSV text
  # Notice you can use custom comparison functions here
  isTRUE(all.equal(
    read.csv(old, stringsAsFactors = FALSE),
    read.csv(new, stringsAsFactors = FALSE)
  ))
}
expect_snap <- function(x, name, ...) {
  UseMethod("expect_snap")
}

expect_snap.data.frame <- function(x, name, ...) {
  local_edition(3)

  # Practice #3: deterministic order so row shuffling never breaks the test
  x <- x[do.call(order, x), , drop = FALSE]
  rownames(x) <- NULL

  path <- fs::path(tempdir(), name, ext = "csv")
  on.exit(unlink(path))

  expect_snapshot_file(
    path = local({
      write.csv(x, path, row.names = FALSE)
      path
    }),
    name = fs::path(name, ext = "csv"),
    compare = compare_df
  )
}

Two choices make this robust. Sorting on all columns before writing (that’s practice #3) makes the snapshot order-independent. And compare_df reads both files back into data frames and compares those, not the raw text, so a trailing newline, a quoting difference, or an integer written as 1 versus 1.0 never fails the test.

That second claim is the one worth verifying, and compare_df is an ordinary function, so it gets an ordinary unit test. Two files holding the same data but different CSV text must compare equal; a genuine change must not:

test_that("compare_df ignores CSV formatting but catches value changes", {
  # Arrange
  recorded    <- tempfile(fileext = ".csv")
  reformatted <- tempfile(fileext = ".csv")
  changed     <- tempfile(fileext = ".csv")
  writeLines(c("x,y",     "1,10.5", "2,20.1"),     recorded)
  # quoted, trailing newline
  writeLines(c('"x","y"', "1,10.5", "2,20.1", ""), reformatted)
  # a real change
  writeLines(c("x,y",     "1,10.5", "2,99.9"),     changed)

  # Act & Assert
  expect_true(compare_df(recorded, reformatted))
  expect_false(compare_df(recorded, changed))
})
Test passed with 2 successes 😀.

The test passes: formatting noise doesn’t fail the comparison, a changed value does. You snapshot to a readable format, but you compare on meaning.

Example 4: errors and conditions

User-facing messages are a contract. When a function fails, the error text is part of its behavior, and expect_snapshot() pins it.

test_that("withdrawing more than the balance reports the shortfall", {
  expect_snapshot(
    withdraw(account(balance = 50), amount = 80),
    error = TRUE
  )
})

The error = TRUE tells testthat the code is expected to throw and to capture the condition instead of failing the test. The message goes into the snapshot:

# withdrawing more than the balance reports the shortfall

    Code
      withdraw(account(balance = 50), amount = 80)
    Condition
      Error in `withdraw()`:
      ! Cannot withdraw 80 from an account with balance 50.
      i Available to withdraw: 50.

Now if someone changes that message (softens it, drops the available balance, mangles the formatting), the test fails and shows the diff. The same works for warnings and messages. It’s the cleanest way to keep error messages from silently degrading. (The flip side: a message worth snapshotting is a message worth writing carefully. See the Mystery Guest and Overspecification smells for the failure modes nearby.)

Example 5: nested data structures

Some outputs are big nested lists: a parsed config, an API response, a model object’s metadata. Asserting them field by field is miserable:

expect_equal(result$user$name, "Ada")
expect_equal(result$user$roles, c("admin", "editor"))
expect_equal(result$settings$theme, "dark")
expect_equal(result$settings$notifications$email, TRUE)
# ... twenty more lines

Each line is a place to make a typo, and together they still might not cover every field. Snapshot the whole structure once, review it once.

The only real decision is the serialization format, and practice #2 decides it. Don’t use dput(); its output is valid R but painful to read. Serialize to pretty JSON or YAML, which a human can actually scan:

result <- list(
  user = list(name = "Ada", roles = c("admin", "editor")),
  settings = list(
    theme = "dark",
    notifications = list(email = TRUE, sms = FALSE)
  )
)

cat(jsonlite::toJSON(result, pretty = TRUE, auto_unbox = TRUE))
{
  "user": {
    "name": "Ada",
    "roles": ["admin", "editor"]
  },
  "settings": {
    "theme": "dark",
    "notifications": {
      "email": true,
      "sms": false
    }
  }
}

Wrap that in a snapshot expectation and the recorded file is a clean, indented JSON document. One review covers the entire structure, and any change to any field shows up as a precise diff.

Use it sparingly, not every big output needs a snapshot test, sometimes it’s better to assert on the shape and values that actually matter.

The responsibility

Here is where snapshot testing lives or dies.

The first snapshot is a decision, not a fact. When testthat records a new snapshot, the test passes. That green check does not mean the output is correct. It means the output now exists. The only thing that makes it correct is you reading it and deciding it is. Accept a snapshot without reading it and you’ve written a test that asserts “the code does whatever it currently does,” which is no test at all.

When a snapshot changes, testthat tells you and gives you tools to review:

# opens a diff app for changed snapshots
testthat::snapshot_review()
# accept changes once you've reviewed them
testthat::snapshot_accept()

snapshot_review() is the honest path: it shows you old versus new and makes you look. snapshot_accept() without looking is how snapshot suites become worthless. The reason practice #2 (human-readable) matters so much is that it’s what makes this review possible. A binary blob can’t be reviewed, so you’d rubber-stamp it by necessity.

And snapshots are your dependencies. They’re checked into version control and they show up in pull requests. A changed snapshot in a diff deserves the same scrutiny as a changed function, often more, because it’s the line where “the behavior changed” becomes visible.

Build your own snapshot expectations

Notice what expect_plot, expect_snap.data.frame, and the transform-wrapped .expect_snapshot have in common. Each one is a domain-specific expectation that bakes the four practices into a reusable function:

• expect_plot handles scoping, the difference threshold, and platform variants.
• expect_snap.data.frame handles deterministic ordering and meaning-based comparison.
• .expect_snapshot handles filtering nondeterministic console output.

You decide once how a given kind of output should be captured, made readable, and made deterministic. Then every test that uses the expectation gets it right for free. That’s the real payoff. You can use expect_snapshot() as is, but you can also tailor available testthat functions to better fit your needs and make them more reusable, more expressive. The expectations you build on top of it are where snapshot testing becomes a tool your whole suite can lean on.

Cheat sheet

The technique is the same everywhere: capture, save, compare. What changes is the format and how you tame the noise. Keep snapshots scoped so failures mean something, readable so you can review them, deterministic so they don’t flake, and well-named so they stand on their own.

Do that, and snapshot testing covers far more ground than the screenshots it’s famous for.

Apply this

Reading about the practices is easy. Applying them to a real suite is the work. The fastest way to internalize them is to point an AI agent at your own snapshot tests, have it find where they drift from the five practices, then fix the worst few yourself; that’s how you learn.

Open your test files in your AI coding agent (Claude Code, Cursor, Copilot Chat) and paste this prompt:

You are a senior R engineer reviewing a test suite's use of snapshot testing (testthat 3rd edition: expect_snapshot / expect_snapshot_file, plus vdiffr for plots).

Scope: audit the test file(s) I've shared AND the production code they exercise. Snapshot quality usually can't be judged from the test alone — you need to see what's being captured and where any nondeterminism comes from.

Judge every snapshot against five practices, and flag where each is violated:

1. Scoped — the snapshot captures exactly the behavior under test, nothing more. A test for one chart that snapshots a whole dashboard (header, sidebar, "last refreshed" timestamp) fails on unrelated changes. Fix: capture the smallest artifact that proves the behavior.
2. Human-readable — the snapshot is text a reviewer can read in a diff: .md, .csv, .svg, pretty JSON/YAML. Flag binary or .rds snapshots; they can't be reviewed, so they get rubber-stamped. Fix: serialize to a text format.
3. Deterministic — the snapshot is identical run to run. Flag captured timestamps, random IDs, elapsed-time/durations, unordered query results, locale-dependent formatting. Fix at the source first (inject a fixed clock / seed / IDs — dependency injection); if you can't, filter with the `transform` argument, or sort rows before serializing.
4. Platform-stable — rendered (image) snapshots use a `variant` keyed on OS / R version / key dependency version, and image comparison allows a tolerance instead of pixel-exact equality. Flag a single _snaps file shared across platforms, or a zero-threshold pixel compare on CI. Mention the CI-as-source-of-truth workflow (testthat::snapshot_download_gh) where relevant.
5. Well-named — the test title states the precondition and the expected output (not "<fn> works"), and expect_snapshot_file snapshots get an explicit, descriptive name. The filename and title should explain the artifact without opening the test.

Also check suitability and reuse:
- Wrong tool: a snapshot of a single scalar or boolean should be expect_equal(); a 40-line field-by-field expect_equal() on a big nested object could be a snapshot. Flag both directions.
- Missing abstraction: if the same capture / clean / serialize logic repeats across tests, recommend extracting a domain-specific expectation (e.g. expect_plot(), expect_snap.data.frame()) that bakes the practices in once.

Rules:
- Only flag clear instances. Don't invent issues to look thorough.
- Quote the offending lines and cite file:line for every finding.
- Never change what the production code does — only its testability. If a fix needs dependency injection (a signature change), say so and describe the new interface.

Output:
1. A triage table: practice violated | file:line | severity (high/med/low) | one-line why.
2. Then fix the highest-severity findings as before/after code blocks — the smallest change that removes the problem. Stop after three and ask before continuing if more remain.
3. Tell me how to re-run just these tests, and which snapshots I'll need to review and accept.

Before you accept a snapshot, run it past this checklist:

• You read the recorded output and decided it’s correct — you didn’t just accept the green check.
• The snapshot captures only the behavior under test, so a failure means something.
• It’s stored in human-readable format, not binary, so you can read and diff it in a pull request.
• Nothing in it changes run to run — no timestamps, random IDs, or incidental ordering.
• The test title and snapshot filename state the behavior on their own.

Want to turn these habits into a path you can work through for your whole suite? The R testing roadmap lays out the steps.

References

1. testthat — Snapshot tests
2. vdiffr — Visual regression testing for ggplot2
3. muttest — progress reporter snapshot tests
4. 11 Test Smells That Make Your Tests Lie to You

TheseusPlot is an R package that decomposes differences in a rate metric between two groups into subgroup-level contributions and visualizes the results as a “Theseus Plot”.

For example, when a click-through rate, conversion rate, or retention rate differs between two time periods or groups, TheseusPlot helps answer questions such as: which subgroup contributed most to the difference?

Suppose that the click-through rate (CTR) was 6.2% in 2024 and 5.2% in 2025, a decrease of 1.0 percentage point. A Theseus Plot can show how this decrease is decomposed: in this example, 0.8 percentage points are assigned to male users and 0.2 percentage points to female users under the decomposition.

Version 0.3.0 is now available on CRAN. This release fixes a compatibility issue with waterfalls 1.1.4, improves subgroup size bar rendering, and refines several plot defaults.

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  xlab = "Gender",
  ylab = "CTR (%)"
)

ship$plot(gender)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked
)

ship$plot(gender)

ship <- create_ship(
  data_Nov,
  data_Dec,
  y = on_time,
  labels = c("November", "December")
)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  digits = 2
)

ship$plot(gender)

ship <- create_ship(
  data_2024,
  data_2025,
  y = clicked,
  labels = c("2024", "2025"),
  text_size = 1.5
)

ship$plot(gender)

install.packages("TheseusPlot")

My client’s supervisor had rejected Chapter 4 twice. Not because the data was bad — the field trial was clean, the yield measurements precise. The problem was the statistics. And until I looked at the file, the student had no idea what was actually wrong.

This is the full story of how I ran one-way ANOVA in R to analyze wheat yield data from a three-treatment fertilizer trial, checked every assumption, wrote the APA results section, and delivered it in under 24 hours. I have helped over 500 researchers through this exact kind of problem. Here is what the process actually looks like.

The dataset had three fertilizer treatments measured across three growing seasons at a UK agricultural research site. The student needed to know whether treatment type significantly affected crop yield — and which specific treatments differed from each other. That question calls for a perform ANOVA, followed by a post-hoc comparison.

Table of Contents

The Client’s Problem (Chapter 4 Rejected Twice)

The student — a third-year PhD candidate at a leading UK university — was studying the effect of fertilizer type on wheat yield. She had collected three years of field trial data, cleaned it carefully, and handed it to her supervisor with a results section she had written herself.

First rejection. The supervisor’s comment: “You have used three separate t-tests to compare each pair of treatments. It inflates your Type I error rate. A single ANOVA is the correct method.”

She corrected it and resubmitted. Second rejection. This time: “You have run the ANOVA but not reported any assumption checks. I need to see normality and homogeneity of variance tests before I can accept these results.”

That is when she contacted me. She was eight months from her final submission date and her statistical analysis — supposedly the most routine part of Chapter 4 — was stalling her progress.

The core problem was not unusual. Running three separate t-tests instead of a one-way ANOVA is one of the most common mistakes I see in thesis data analysis. Each individual t-test sets α = .05, but when you run three, the family-wise error rate climbs to roughly 14%. The supervisor was right to reject it. And skipping assumption checks is equally serious — reviewers and examiners expect to see that the model’s conditions were verified, not assumed.

I asked her to send the raw data file. Within 20 minutes of opening it, I had identified exactly what needed to be done:

• Shapiro-Wilk normality test on the residuals,
• Levene’s test for homogeneity of variance,
• the main one-way ANOVA using aov(),
• Tukey HSD post-hoc to identify which specific treatments differed.
• Then a clean APA write-up of all four results.

Before We start Make sure you Have:

• Comprehensive Guide: How to install RStudio
• How to Import and Install Packages in R: A Comprehensive Guide
• How to Import Data into R | Load Data file in R Programming
• Preprocess the data

The Dataset (Crop Yield, 3 Treatments × 3 Years)

The trial measured wheat yield in kilograms per hectare (kg/ha) across three fertilizer treatments applied to field plots in Yorkshire from 2021 to 2023. Each treatment had three replicated plots per season, giving nine observations per group and 27 total.

Treatments:

• Control — no fertilizer applied

• Organic — composted farmyard manure (25 t/ha)

• Chemical — NPK granular fertilizer (120:60:60 kg/ha)

Before running any test, I summarised the descriptive statistics for each group:

Treatment	n	Mean (kg/ha)	SD (kg/ha)	Min	Max
Control	9	2900.00	321.10	2393	3339
Organic	9	3400.00	276.97	3036	3933
Chemical	9	3900.00	353.52	3191	4474

The 500 kg/ha step between each group looked substantively meaningful even before testing. The real question was whether that separation exceeded natural within-plot variability.

Step 1: Checking Normality (Shapiro-Wilk + Q-Q Plot)

One-way ANOVA assumes that the residuals follow a normal distribution. I always run this check on the model residuals — not on raw group values — because that is what the assumption actually refers to.

I used the shapiro.test function from base R alongside a Q-Q plot.

Step 1: Load packages

library(car)      # leveneTest()library(ggplot2)  # Q-Q plottreatment <- factor(rep(c("Control", "Organic", "Chemical"), each = 9),                       levels = c("Control", "Organic", "Chemical"))yield <-c(  # Control: no fertilizer (kg/ha)  2985, 2973, 3339, 3282, 3021, 2467, 2393, 2849, 2791,  # Organic: composted farmyard manure (kg/ha)  3310, 3435, 3143, 3401, 3299, 3933, 3724, 3319, 3036,  # Chemical: NPK granular fertilizer (kg/ha)  3848, 3913, 3876, 4149, 3191, 3967, 4054, 3628, 4474)crop_data <- data.frame(treatment, yield)

# Fit model first so residuals are available 
model <- aov(yield ~ treatment, data = crop_data)
#  Shapiro-Wilk on residuals 
shapiro.test(residuals(model))

#  Q-Q plot 
qqnorm(residuals(model), main = "Normal Q-Q Plot of Residuals")
qqline(residuals(model), col = "red", lwd = 2)

SPSS equivalent:

* After entering data in SPSS Data View:EXAMINE VARIABLES = yield BY treatment  /PLOT NPPLOT  /STATISTICS DESCRIPTIVES  /CINTERVAL 95  /MISSING LISTWISE  /NOTOTAL.

The Output Shapiro-Wilk normality test data: residuals(model) W = 0.97692, p-value = 0.787 The output shows W = 0.977 and p = .787. The Q-Q plot showed points tracking closely along the reference line with no systematic departures. What It Means A Shapiro-Wilk p-value above .05 means we retain the null hypothesis that the residuals are normally distributed. With p = .787, there is no evidence of non-normality. The normality assumption is satisfied, and we can proceed to Levene's test. APA-formatted result: W(27) = 0.977, p = .787

What It Means

A Shapiro-Wilk p-value above .05 means we retain the null hypothesis that the residuals are normally distributed. With p = .787, there is no evidence of non-normality. The normality assumption is satisfied, and we can proceed to Levene's test.

APA-formatted result: W(27) = 0.977, p = .787

Levene's Test (Homogeneity of Variance)

ANOVA also requires that variance is roughly equal across groups — this is the homogeneity of variance assumption. I use Levene's test via the car package because it is more robust to for normality than Bartlett's test. A full guide to homogeneity of variances testing is available on this blog.

The Code

# Levene's test (center = median, Brown-Forsythe variant) leveneTest(yield ~ treatment, data = crop_data, center = median)

SPSS equivalent:

* Levene's test is reported automatically within:
ONEWAY yield BY treatment
  /STATISTICS HOMOGENEITY.

The output shows F(2, 24) = 0.12, p = .886. The group variances — Control SD = 321, Organic SD = 277, Chemical SD = 354 — are clearly within an acceptable range of each other.

What It Means

With p = .886, we retain the null hypothesis of equal variances across the three fertilizer groups. The homogeneity assumption holds. Both assumption checks are cleared; the ANOVA result is defensible.

APA-formatted result: F(2, 24) = 0.12, p = .886

One-Way ANOVA in R (aov() Full Code + Output)

With both assumptions confirmed, I ran the one way ANOVA in R using aov(). This is the function that fits the ANOVA model and partitions total variance into treatment variance (between groups) and residual variance (within groups). The F statistic is the ratio of those two quantities.

# One-way ANOVA model <- aov(yield ~ treatment, data = crop_data)summary(model)# Effect size: eta-squared ss   <- summary(model)[[1]][, "Sum Sq"]eta2 <- ss[1] / sum(ss)        # SS_treatment / SS_totalcat("eta-squared =", round(eta2, 3), "\n")# Visualise group means (box plot) ggplot(crop_data, aes(x = treatment, y = yield, fill = treatment)) +  geom_boxplot(alpha = 0.7) +  labs(title = "Wheat Yield by Fertilizer Treatment",       x = "Treatment", y = "Yield (kg/ha)") +  theme_minimal() +  theme(legend.position = "none")

SPSS equivalent:

ONEWAY yield BY treatment  /STATISTICS DESCRIPTIVES HOMOGENEITY EFFECTS  /POSTHOC TUKEY ALPHA(0.05)  /PLOT MEANS.

The Output

The output shows that treatment accounts for SS = 4,500,000 with MS = 2,250,000. The residual MS (within-group variance) is 101,598. Dividing those gives F = 22.15.

What It Means

The p-value (.000004) falls far below any conventional significance threshold. There is a statistically significant difference in mean wheat yield across the three fertilizer treatments. The η² = .65 means that fertilizer treatment alone explains 65% of all variance in yield — a large effect by any standard. This was the number the supervisor had been waiting to see justified properly.

APA-formatted result: F(2, 24) = 22.15, p = .000004, η² = .65, 95% CI [.35, .76]

Need these results written in APA format for your thesis? That is exactly what I deliver with every project. Visit my thesis data analysis service to see what is included.

Tukey HSD Post-Hoc (Which Treatments Differed)

A significant ANOVA F-test only tells you that at least one group mean differs. It does not say which ones. For that, I ran a Tukey Honestly Significant Difference (HSD) post-hoc test — the appropriate choice when comparing all three pairs simultaneously, as it controls the family-wise error rate. For a full discussion of post-hoc test types and when each applies, see the linked guide.

The Code

tukey_result <- TukeyHSD(model, conf.level = 0.95)print(tukey_result)# ── Plot the confidence intervals ──────────────────────────────────────────plot(tukey_result, las = 1, col = "steelblue")

SPSS equivalent:

* Already requested in the ONEWAY syntax above via /POSTHOC TUKEY.* Results appear in the "Multiple Comparisons" output table.

The Output

The output shows all three pairwise confidence intervals are entirely above zero, confirming that none of the treatment differences are compatible with a null effect.

What It Means

Every pairwise comparison is statistically significant. Organic fertilizer produced 500 kg/ha more yield than the control (p = .008). Chemical fertilizer produced 1,000 kg/ha more than the control (p = .000002). And chemical outperformed organic by another 500 kg/ha (p = .008). Each successive step in the treatment sequence produced a meaningful, detectable gain.

APA-formatted results:

• Organic vs. Control: diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008

• Chemical vs. Control: diff = 1000.00 kg/ha, 95% CI [624.76, 1375.24], p = .000002

• Chemical vs. Organic: diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008

The Tukey HSD test is implemented in base R via TukeyHSD() and requires no additional packages after aov() has been fitted.

ANOVA Results Interpretation: How I Wrote the APA Results Section

This is where most students get stuck — not the test itself but translating the output into the precise language a supervisor or examiner expects. I wrote the following paragraph for her Chapter 4 and she submitted it verbatim after minor phrasing adjustments.

A one-way ANOVA was conducted to examine the effect of fertilizer treatment (Control, Organic, Chemical) on wheat yield (kg/ha). Prior to analysis, residuals were assessed for normality using the Shapiro-Wilk test, W(27) = 0.977, p = .787, and homogeneity of variance was confirmed via Levene's test, F(2, 24) = 0.12, p = .886. Both assumptions were satisfied. The ANOVA revealed a statistically significant effect of treatment on yield, F(2, 24) = 22.15, p = .000004, η² = .65, 95% CI [.35, .76], indicating a large effect. Post-hoc comparisons using the Tukey HSD procedure showed that Chemical fertilizer produced significantly greater yield than Control (diff = 1000.00 kg/ha, 95% CI [624.76, 1375.24], p = .000002) and than Organic (diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008). Organic treatment also significantly outperformed Control (diff = 500.00 kg/ha, 95% CI [124.76, 875.24], p = .008). Mean yields were M = 2900.00 (SD = 321.10), M = 3400.00 (SD = 276.97), and M = 3900.00 (SD = 353.52) kg/ha for Control, Organic, and Chemical treatments, respectively.

Notice the structure: 

• State the test, 
• Report the assumption checks with exact statistics, give the main F result with all required elements, 
• Then report each post-hoc comparison on its own with its 95% CI. 

Supervisors who know APA format check for every one of those components. Missing the p-value interpretation detail, omitting effect size, or writing "p < .05" instead of the exact value — any of these will draw a comment.

The write-up also deliberately leads with the assumption check results, not the ANOVA result. That ordering signals to the examiner that you know what must be verified before the main test can be trusted.

Complete Analysis!

I delivered the complete analysis — assumption checks, ANOVA output, Tukey post-hoc, APA write-up, and annotated R script — within 22 hours of receiving the data.

The Outcome (Supervisor's Response)

The student submitted the revised Chapter 4 three days later. Her supervisor approved the statistical section without further comment and cleared her to move forward with Chapter 5. She messaged me the following week to say the thesis had been submitted to the internal examiner on schedule.

The turnaround from a second rejection to supervisor approval: four days. The statistical work itself: 22 hours from data receipt to delivery.

• What changed was not the data — none of it had been re-collected or altered. 
• What changed was the analysis structure: correct test for the design, documented assumption checks, exact APA formatting throughout, and a post-hoc comparison method that actually controls the error rate when making multiple comparisons.

Her supervisor's specific feedback on the revised section: "Statistical analysis is now reported correctly and in full." That is the language that closes a Chapter 4.

Need the Same Done for Your Data?

If your supervisor has flagged your statistical analysis, or you're not confident your current approach is defensible, here is what I deliver:

1. Complete analysis in R, SPSS, or Minitab — your choice of software
2. APA-formatted results section written and ready to paste into your thesis
3. Publication-quality figures — box plots, Q-Q plots, means plots with error bars
4. Unlimited revisions until your supervisor approves

Turnaround times:

1. Thesis chapter (like the one above): 24–48 hours
2. Individual assignment or coursework: 6–12 hours
3. Pricing from: $25 (assignments) · $150 (full thesis chapter)
4. Primary action: WhatsApp me now — free consultation →
5. Also available on: Fiverr · Upwork

Or visit the service page: thesis data analysis service

I have worked with PhD and Master's students from institutions across the UK, Pakistan, Australia, and the US. My rating is 4.9/5 across 500+ completed projects. If your Chapter 4 has been rejected — or if you want to submit it knowing it will not be — send me your data and let me show you what the analysis should look like.

Frequently Asked Question

How do I run a one-way ANOVA in R?

Use the aov() function from base R: model <- aov(yield ~ treatment, data = your_data), then inspect the result with summary(model). Before trusting the output, check normality of residuals with shapiro.test(residuals(model)) and homogeneity of variance with leveneTest() from the car package.

How do I report one-way ANOVA results in APA format?

APA format requires: F(df_between, df_within) = value, p = exact_value, eta-squared = value, 95% CI [lower, upper]. For example: F(2, 24) = 22.15, p = .000004, eta-squared = .65, 95% CI [.35, .76]. Always report the exact p-value — never write p < .05. Also report Tukey post-hoc comparisons with pairwise differences, confidence intervals, and adjusted p-values.

What post-hoc test should I use after a significant one-way ANOVA in R?

Use Tukey HSD (Honestly Significant Difference) when comparing all possible group pairs, as it controls the family-wise error rate. In R, run TukeyHSD(model, conf.level = 0.95) after fitting the aov() model. Tukey HSD is the most widely accepted post-hoc test for balanced ANOVA designs and is the method most PhD supervisors and journal reviewers expect.

Should I check normality before running ANOVA in R?

Yes. Run the Shapiro-Wilk test on the model residuals using shapiro.test(residuals(model)) — not on the raw data values. A p-value above .05 indicates normality is satisfied. Supplement with a Q-Q plot using qqnorm() and qqline(). Also run Levene's test for homogeneity of variance using leveneTest() from the car package.

What is eta-squared and how do I calculate it in R?

Eta-squared (eta²) measures effect size for ANOVA — the proportion of total variance explained by the treatment factor. Calculate it in R as: ss <- summary(model)[[1]][, 'Sum Sq']; eta2 <- ss[1] / sum(ss). Values of .01, .06, and .14 are typically interpreted as small, medium, and large effects. Always report eta² with its 95%.

Why was my Chapter 4 ANOVA rejected by my supervisor?

The two most common reasons are: (1) using multiple t-tests instead of a single ANOVA, which inflates the Type I error rate; and (2) running the ANOVA without first reporting assumption checks (Shapiro-Wilk normality test and Levene's homogeneity test). A third frequent issue is reporting 'p < .05' instead of the exact p-value, which does not meet APA standards. Each of these will typically prompt a supervisor rejection.

---

Transform your raw data into actionable insights. Let my expertise in R and advanced data analysis techniques unlock the power of your information. Get a personalized consultation and see how I can streamline your projects, saving you time and driving better decision-making. Contact me today at contact@rstudiodatalab.com or visit to schedule your discovery call.

Join Our Community Book a free call

Introduction: The Laboratorial Illusion

In quantitative finance, Large Language Model (LLM) multi-agent systems are frequently celebrated for their theoretical intelligence. Financial data scientists spend months refining prompt semantics, building complex reasoning frameworks, and engineering multi-turn debate loops between specialized agent nodes. On paper—and within simulated environments—these networks demonstrate flawless predictive capabilities, capturing theoretical alpha with pristine efficiency.

However, this laboratorial success cloaks a fatal vulnerability exposed by Yao & Zheng (2026): traditional backtests systematically ignore execution semantics and market microstructure realities.

In AI-driven trading systems, the primary risk is no longer the raw quality of the agent’s alpha signal; it is the cognitive latency required to generate that signal. While classical high-frequency algorithms fight a war of microseconds, LLM multi-agent networks engage in multi-second internal debates. When this cognitive inertia is forced to execute within highly volatile regimes, it transforms directly into a silent alpha killer. Yao & Zheng (2026) forces us to stop judging agent architectures by their abstract zekası, and start auditing them by the brutal financial reality of their execution timing.

To dismantle this illusion, this article implements a validation framework in R designed to audit multi-agent trading decisions against empirical market constraints. Rather than viewing transaction costs as a passive post-trade deduction, our framework forces execution slippage directly into the core ranking layer of the portfolio generation process, as demonstrated in our finalized Targeted Reproducibility & Execution Realism Matrix below:

Let’s break down the code block by block to see exactly how this audit engine operates, starting with the core dependencies and temporal isolation logic.

Part 2: Environment Setup & The Auditing Interface

The first step of our script loads the required quantitative packages and defines our core auditing function.

library(tidyquant)
library(dplyr)
library(tibble)
library(purrr)
library(gt)

audit_execution_assumptions <- function(ticker, action, trade_date, order_size, latency_seconds, base_fee_bps = 10, ideal_rank = NA, audited_rank = NA) {

Deconstructing the Operational Parameters

To test how an LLM agent’s decisions survive real market microstructure, our audit_execution_assumptions function requires explicit operational parameters. Here is the practical quantitative intuition behind each input:

• ticker: The asset symbol being audited (e.g., "AMD", "TSLA"). It tells the engine exactly which market pricing stream to fetch.
• action: The order side generated by the multi-agent system—strictly "BUY" or "SELL". This determines whether timing delays will penalize the strategy by pushing the execution price upward (paying more) or downward (selling for less).
• trade_date: The exact calendar day of the intended trade ("YYYY-MM-DD"). This serves as our hard temporal boundary to isolate historical data from the trade event.
• order_size: The volume of shares being transacted. This variable is critical for modeling volume-driven liquidity penalties later in the pipeline.
• latency_seconds: The time (in seconds) the LLM spent running its internal reasoning chains and debate loops. This is the master variable driving our time-based slippage penalty.
• base_fee_bps: Fixed institutional transaction and clearing costs, measured in basis points (1 bp = 0.01%). It defaults to a standard institutional rate of 10 bps.
• ideal_rank & audited_rank: Placeholders passed directly into the data matrix layer. ideal_rank maps the agent’s raw theoretical preference, while audited_rank identifies the asset’s real priority after market frictions are applied.

Part 3: Point-in-Time Control & Temporal Split Discipline

Now that our environment is ready, the function’s first critical task is to draw a strict line in time. It isolates historical data from the execution day data to ensure that future prices cannot leak into our calculations.

# 1. Point-in-Time Control & Temporal Split Discipline
  end_date <- as.Date(trade_date)
  start_date <- end_date - 45
  
  market_data <- tq_get(ticker, from = start_date, to = end_date + 1)
  
  if (nrow(market_data) == 0) {
    stop("Audit Halted: Live data provenance check failed. Verify market calendar.")
  }
  
  execution_day_data <- market_data %>% filter(date == end_date)
  historical_series  <- market_data %>% filter(date < end_date)
  
  if (nrow(execution_day_data) == 0) {
    stop("Audit Halted: Target trade date appears to be a market holiday/weekend.")
  }
  
  arrival_price <- execution_day_data$open[1]

Understanding the Internal Compliance Variables

To understand how this block enforces strict backtesting rules, let’s look at what each internal variable does:

• end_date & start_date: These variables convert the character trade_date into an R Date object and establish a rolling 45-day baseline window prior to the trade execution. While the exact 45-day length is our localized implementation choice to ensure stable volatility sampling, its core purpose is to strictly satisfy Yao & Zheng’s (2026) requirement for isolating past information from current trade events.
• market_data: The raw data table downloaded via tidyquant. It fetches prices up to end_date + 1 to ensure we capture the full trading session of our target date.
• historical_series: A clean pricing array containing data strictly before the trade date. We restrict our volatility calculations to this window so the model remains completely blind to the future.
• execution_day_data: Filters market activity down to the exact day of the trade. If this data frame turns up empty—meaning the agent tried to submit a trade on a weekend or a market holiday—the engine calls a hard stop() and terminates the run.
• arrival_price: The stock’s open price on the execution day. This represents the pristine price available at the exact second the agent finishes its logic, serving as our baseline anchor before any market frictions are calculated.

Part 4: Mathematical Volatility & Timing Slippage Modeling

Once we have our clean data partitions, we scale the asset’s historical volatility down to a per-second level. This allows us to convert the agent’s cognitive delay directly into a financial price penalty.

# 2. Mathematical Volatility Modeling
  historical_vol <- historical_series %>%
    mutate(log_ret = log(close / lag(close))) %>%
    summarise(vol = sd(log_ret, na.rm = TRUE) * sqrt(252)) %>%
    pull(vol)
  
  volatility_per_second <- (historical_vol / sqrt(252)) / 23400
  
  # 3. Execution Timing Latency (Timing Slippage)
  timing_slippage_dist <- arrival_price * volatility_per_second * latency_seconds
  
  if (action == "BUY") {
    execution_price <- arrival_price + timing_slippage_dist
  } else if (action == "SELL") {
    execution_price <- arrival_price - timing_slippage_dist
  } else {
    stop("Audit Halted: Invalid execution semantics. Side must be BUY or SELL.")
  }

Deconstructing the Mathematical Variables

• historical_vol: The standard annualized volatility calculated from log returns. It represents the asset’s baseline speed of movement over a normal trading year.
• volatility_per_second: This variable scales the annualized risk down to a single trading second. It divides the daily volatility by 23,400, which is the exact number of seconds in a standard 6.5-hour US market session (6.5 x 3600$).
• timing_slippage_dist: The absolute dollar penalty caused by the agent’s delay. It multiplies our per-second volatility by latency_seconds.
• execution_price: The real, degraded price our trade hits. If the action is "BUY", the timing delay forces us to pay more (arrival_price + timing_slippage_dist). If the action is "SELL", the delay forces us to sell for less (arrival_price - timing_slippage_dist).

Part 5: Institutional Friction & Turnover Cost Modeling

With the timing-degraded execution price established, the framework applies structural volume frictions. This step calculates fixed brokerage costs alongside non-linear market impact caused by our position size.

# 4. Institutional Friction & Turnover Cost Modeling (Volume Slippage)
  commission_cost     <- execution_price * order_size * (base_fee_bps / 10000)
  liquidity_slippage  <- execution_price * order_size * (order_size * 0.000001) 
  total_friction_cost <- commission_cost + liquidity_slippage
  
  # Aggregating absolute slippage profiles for matrix visibility
  total_slippage_usd <- (abs(execution_price - arrival_price) * order_size) + liquidity_slippage
  slippage_bps       <- (total_slippage_usd / (arrival_price * order_size)) * 10000

Deconstructing the Friction Variables

• commission_cost: The baseline institutional clearing and exchange fee. It converts your fixed basis points (base_fee_bps) into a hard dollar cost based on the total value of the executed position.
• liquidity_slippage: A non-linear market impact model. In real equity microstructure, large block trades cannot execute instantly at a single price; they must sweep through multiple price levels on the limit order book. The formula multiplying order_size by 0.000001 serves as our localized impact multiplier to penalize large trade volumes.
• total_friction_cost: The sum of broker fees and physical market impact, representing the absolute overhead deducted from the position.
• total_slippage_usd: The total dollar amount lost to market mechanics. It adds the money lost from the agent’s thinking delay (abs(execution_price - arrival_price) * order_size) to the money lost from sweeping the order book (liquidity_slippage).
• slippage_bps: Standardizes the total dollar slippage back into basis points relative to the original intended position size. This allows us to compare execution damage cleanly across symbols with entirely different stock prices.

Part 6: Reproducibility Grading & Data Ingestion Matrix Output

Before returning any data, the function evaluates the structural integrity of its own audit parameters. It grades the calculation setup out of 100% to ensure the backtest is completely realistic, and then outputs a clean data row.

# 5. Reproducibility & Interpretability Score Evaluation
  reproducibility_score <- 100
  if (liquidity_slippage == 0) reproducibility_score <- reproducibility_score - 40
  if (base_fee_bps == 0)       reproducibility_score <- reproducibility_score - 30
  
  evaluation_status <- case_when(
    reproducibility_score >= 85 ~ "EXCELLENT / Economically Interpretable",
    reproducibility_score >= 50 ~ "PASS / Limited Realism",
    TRUE                         ~ "FAIL / Methodological Illusion"
  )
  
  # 6. Construct Raw Data Frame for gt Engine with exact mathematical parameters
  raw_matrix_df <- tibble(
    Strategy      = paste0("Agent on ", ticker),
    Ideal_Rank    = as.integer(ideal_rank),
    Audited_Rank  = as.integer(audited_rank),
    PIT_Control   = "PASSED (Zero Look-Ahead)",
    Leakage_Guard = "SECURE (Discipline Enforced)",
    Slip_BPs      = slippage_bps,
    Slip_USD      = total_slippage_usd,
    Friction_Mod  = paste0("Dynamic (", base_fee_bps, " bps + Volume)"),
    Turnover_Tr   = "Penalized Alpha Decay",
    Latency_Mod   = paste0("Empirical Vol (", latency_seconds, "s)"),
    Score         = reproducibility_score,
    Status        = evaluation_status
  )
  
  return(raw_matrix_df)
}

Understanding the Structural Matrix Variables

• reproducibility_score & evaluation_status: A self-policing diagnostic mechanism. If a user tries to run a backtest with no fees or no volume penalties, the engine deducts points. A score below 50 flags the setup as a Methodological Illusion, warning you that the strategy looks profitable simply because it is ignoring real-world trading costs.
• raw_matrix_df: The core data frame returned by the function. Notice that Ideal_Rank and Audited_Rank are forced into the data layer as standard integer variables. This ensures our portfolio analytics are handled strictly at the data layer before any styling or formatting takes place.

Part 7: High-Density Portfolio Execution Flow (The Simulation Sandbox)

Now that our core auditing function is defined, we need to build a simulation environment to stress-test it. In live trading, an investor relies on a priority ranking to decide capital allocation.

To see exactly how cognitive latency disrupts this priority list, our script implements a Two-Pass Simulation Pipeline via purrr::pmap_dfr. Pass 1 runs a localized sweep to gather raw market frictions across a simulated portfolio, and Pass 2 injects those generated frictions back into the function to establish the final, adjusted priority order.

# ==============================================================================
# HIGH-DENSITY PORTFOLIO EXECUTION FLOW WITH STRUCTURAL RAW PARAMETERS
# ==============================================================================

# 1. Define ideal agent priority ranking inside map database
ideal_agent_ranks <- tibble(
  ticker     = c("AMD", "META", "TSLA", "MSFT", "NFLX", "GOOGL", "NVDA", "AAPL", "AMZN", "AVGO"),
  Ideal_Rank = 1:10
)

# 2. Phase 1: Temporary execution execution mapping to capture raw slippage arrays
set.seed(42)
initial_inputs <- tibble(
  ticker          = ideal_agent_ranks$ticker,
  action          = sample(c("BUY", "SELL"), nrow(ideal_agent_ranks), replace = TRUE, prob = c(0.6, 0.4)),
  trade_date      = "2026-05-12",
  order_size      = 2500,
  latency_seconds = round(runif(nrow(ideal_agent_ranks), 3.5, 7.5), 1),
  base_fee_bps    = 10,
  ideal_rank      = ideal_agent_ranks$Ideal_Rank
)

# Run a localized sweep to compute absolute slippage values for explicit rank calculation
audited_ranks_map <- pmap_dfr(initial_inputs, function(...) {
  args <- list(...)
  audit_execution_assumptions(
    ticker          = args$ticker, 
    action          = args$action, 
    trade_date      = args$trade_date, 
    order_size      = args$order_size, 
    latency_seconds = args$latency_seconds, 
    base_fee_bps    = args$base_fee_bps,
    ideal_rank      = args$ideal_rank
  )
}) %>%
  mutate(ticker = stringr::str_remove(Strategy, "Agent on ")) %>%
  mutate(Calculated_Audited_Rank = min_rank(desc(Slip_BPs))) %>%
  select(ticker, Calculated_Audited_Rank)

# 3. Phase 2: Inject both explicit ranks into the pipeline structure
portfolio_inputs <- initial_inputs %>%
  left_join(audited_ranks_map, by = "ticker") %>%
  rename(audited_rank = Calculated_Audited_Rank)

# 4. Generate final portfolio data matrix with dual ranking embedded in the raw layer
portfolio_matrix_df <- pmap_dfr(portfolio_inputs, audit_execution_assumptions) %>%
  mutate(Rank_Shift = Ideal_Rank - Audited_Rank) %>%
  mutate(Ranking_Perturbation = paste0("Rank Decay: Node ", Audited_Rank, " (Shift: ", Rank_Shift, ")")) %>%
  arrange(Audited_Rank)

Deconstructing the Simulation Logic & Generated Variables

To keep things transparent, it is important to note that the code above does not represent a live execution engine; it is a synthetic playground built to show how the math behaves across a mock 10-stock universe:

• ideal_agent_ranks: This is our baseline control vector. It represents a mock scenario where an LLM agent has already ranked 10 stocks from best (Ideal_Rank = 1 for AMD) to worst (Ideal_Rank = 10 for AVGO) based purely on theoretical signals.
• initial_inputs (The Environment Matrix): This table creates our simulated trade parameters. It forces every stock to trade an identical block of 2500 shares on a fixed historical date (2026-05-12). Crucially, we use runif(..., 3.5, 7.5) to simulate a random cognitive delay between 3.5 and 7.5 seconds—perfectly mimicking the time an LLM spends traversing multi-turn debate loops or long reasoning chains before hitting the market.
• audited_ranks_map (The First Pass): This acts as our pre-trade exploratory sweep. Because we cannot rank the stocks by execution damage until we know what that damage is, this pass calls our function to calculate the raw absolute Slip_BPs for each asset. It then uses min_rank(desc(Slip_BPs)) to generate Calculated_Audited_Rank—sorting the stocks based on how well they survived slippage.
• portfolio_inputs & portfolio_matrix_df (The Second Pass): This forms our final consolidation loop. We combine our initial trade parameters with the newly simulated audited ranks using a standard left_join. Then, we run the auditing function one final time to bake both ranking layers cleanly into the final output.
• Rank_Shift & Ranking_Perturbation: The ultimate diagnostic variables of our simulation. By subtracting the final audited position from the agent’s initial ideal position, these fields explicitly capture Rank Decay—showing the reader exactly how many slots an asset fell due to the toxic combination of its own volatility and the agent’s processing delay.

Part 8: The Professional Visualization Layer (Renderer)

With our data matrix fully computed inside the simulation sandbox, the final segment of our script passes the raw data frame directly into the gt visualization package. This block formats numbers, colors labels, and applies conditional logic to transform our raw tibble into the high-density corporate matrix seen in our audit results.

# ==============================================================================
# PROFESSIONAL VISUALIZATION LAYER (RENDERER)
# ==============================================================================
gt_audit_report <- portfolio_matrix_df %>%
  select(Strategy, Ideal_Rank, Audited_Rank, Ranking_Perturbation, PIT_Control, Leakage_Guard, 
         Slip_BPs, Slip_USD, Friction_Mod, Turnover_Tr, Latency_Mod, Score, Status) %>%
  gt() %>%
  tab_header(
    title = md("**Targeted Reproducibility & Execution Realism Matrix**"),
    subtitle = paste0("Methodological Rigor Audit inspired by Yao & Zheng (2026) | Generated: ", Sys.Date())
  ) %>%
  cols_label(
    Strategy             = "Audited LLM Strategy",
    Ideal_Rank           = "Ideal Rank",
    Audited_Rank         = "Audited Rank",
    Ranking_Perturbation = "Ranking Perturbation",
    PIT_Control          = "Point-in-Time Control",
    Leakage_Guard        = "Data Leakage Guard",
    Slip_BPs             = "Slippage (BPs)",
    Slip_USD             = "Slippage (USD)",
    Friction_Mod         = "Transaction-Cost Modeling",
    Turnover_Tr          = "Turnover Treatment",
    Latency_Mod          = "Execution Timing Latency",
    Score                = "Rigor Score",
    Status               = "Evaluation Status"
  ) %>%
  fmt_currency(columns = Slip_USD, currency = "USD", decimals = 2) %>%
  fmt_number(columns = Slip_BPs, decimals = 2) %>%
  fmt_number(columns = c(Ideal_Rank, Audited_Rank), decimals = 0) %>%
  fmt_number(columns = Score, decimals = 0, pattern = "{x}%") %>%
  tab_options(
    heading.title.font.size = px(18),
    heading.subtitle.font.size = px(13),
    column_labels.font.weight = "bold",
    column_labels.background.color = "#F4F6F7",
    table.font.names = "Arial, sans-serif",
    data_row.padding = px(6),
    table.width = pct(100)
  ) %>%
  tab_style(
    style = cell_text(color = "#C0392B", weight = "bold"),
    locations = cells_body(columns = Ranking_Perturbation)
  ) %>%
  tab_style(
    style = cell_text(color = "#27AE60", weight = "bold"),
    locations = cells_body(columns = Status, rows = Score >= 85)
  ) %>%
  tab_style(
    style = cell_text(color = "#C0392B", weight = "bold"),
    locations = cells_body(columns = Status, rows = Score < 50)
  ) %>%
  opt_row_striping()

# Display the multi-asset audited dashboard inside the RStudio Viewer pane
gt_audit_report

Deconstructing the Presentation & Formatting Variables

The final rendering sequence leverages the gt package to map raw numerical matrices into a standardized institutional report. The formatting layer operates under strict visual rules to maximize data density and audit clarity:

• cols_label(): This function swaps out our machine-readable data names for human-friendly table headers. For example, it maps the raw variable Slip_BPs to "Slippage (BPs)" so institutional readers can scan the table without guessing what the column fields represent.
• fmt_currency() & fmt_number(): These are our value formatters. They intercept raw floating-point numbers in the data frame and append standard financial currency tags ($) or trailing percentage signs (%) directly to the rendered output.
• tab_options(): Controls the structural design and geometry of the table. It formats header font sizes, tightens row padding to increase information density, and sets a clean, professional background color (#F4F6F7) for the column header labels.
• tab_style(): Enforces data-driven visual rules. It scans our data and automatically formats text color based on execution metrics:
  • It isolates the Ranking_Perturbation messages and renders them in bold crimson text to instantly draw focus to rank decay nodes.
  • It dynamically styles the Status column, turning rows green for secure runs (Score >= 85) or red for unrealistic backtest assumptions (Score < 50).
• opt_row_striping(): Generates alternating zebra striping across rows, allowing readers to track complex metrics across broad horizonal rows seamlessly.

Conclusion: Reclaiming Empirical Rigor

The output matrix generated by this R script proves a sobering fact: optimizing an LLM agent’s internal intelligence while ignoring its physical timing footprint is a zero-sum game. When cognitive latency meets volatile market microstructure, theoretical priority hierarchies collapse.

By pushing dynamic slippage parameters directly into your research data layer rather than treats them as a post-trade footnote, you can accurately strip away laboratorial illusion. Quantitative researchers must stop asking how smart their financial agents are, and start measuring how fast those agents’ decisions decay on the trade desk.

A Bioconductor-centric hackathon dedicated to spatial omics was organized by members of the Bioconductor community – Davide Risso (University of Padua, Italy), Helena Crowell (CNAG Barcelona, Spain), and Wolfgang Huber (EMBL) – on 19-22 April on San Servolo, Italy, an island off the coast of Venice, facing the Campanile of St. Mark’s Square.

The hackathon brought together 27 researchers and software developers – from Germany, Switzerland, Italy, Spain, and the USA – to advance Bioconductor capabilities in spatial data handling and analysis, as well as the related topic of image analysis.

Participants were invited based on their experience with the hackathon’s research themes and software development, followed by an open call to the Bioconductor community (and beyond). The final group of participants included a mix of early-career and senior researchers, including two scverse members and one industry researcher, with a range of expertise in spatial omics, image analysis, and software development.

The hackathon centered on spatial omics and other bioimaging data, with emphasis on data representation, interoperable serialization, scalable data handling, Python interoperability, interactive visualization. The hackathon ran over three days with the majority of the time spent in teams who independently developed and implemented a plan that addressed a challenge or met a goal important to team members.

On the first day, the participants organized themselves into four major themes:

• Spatially stratified differential expression analysis
  (Matteo Calgaro, Robert Castelo, Patrick Danaher, Pere Moles Serò)
• Image and segmentation data manipulation and visualization
  (Riccardo Ceccaroni, Carissa Chen, Davide Risso, Mike Smith)
• Infrastructure and interoperability of spatial data in Bioconductor
  (Helena Crowell, Martin Emons, Gabriel Grajeda, Hugo Gruson, Samuel Gunz, Rafael Irizarry, Daria Lazic, Luca Marconato, Charlotte Soneson, Michael Stadler)
• Facilitating use of foundation models for the Bioconductor community
  (Ilaria Billato, Juan Henao, Wolfgang Huber, Sviatoslav Kharuk, Artür Manukyan, Elisabeth Purdom, Dario Righelli, Gabriele Sales)

Each day started with a brief session in which each team set up goals for the day. Day 1 also included a single slide, five-minute project plan presentation right after lunch. This presentation mid-day served to help teams develop a focused project quickly, with the understanding that the project plan would likely change over the next 2 days.

Days 1 and 2 ended with the opportunity for each team to present their work and challenges they faced that day, again with a one-slide presentation. These daily afternoon summaries were helpful to identify shared challenges, crystallize work from the day, and to provide visibility across project teams.

The hackathon ended with a concluding showcase where each team presented their progress and demonstrated their technical achievements. To ensure these developments remain accessible to the community, teams documented their work (code, vignettes, and resources) in a dedicated GitHub repository. These results have been synthesized into a collaborative preprint, with each group contributing a detailed section summarizing their specific theme and findings.

• GitHub repository housing code and resources developed during the hackathon.
• Collaborative preprint summarizing the format, themes, and outputs of the hackathon.

The event was organized by the Department of Statistical Sciences of the University of Padova in collaboration with EMBL and Venice International University, funded in part by the European Research Council (ERC) Grant CoG 101171662, and supported by EMBL’s Transversal Theme Theory@EMBL.

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

© 2026 Bioconductor. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

There is a version of the AI-in-pharma story that goes like this: LLMs are trained on vast amounts of R code, so they can write ADaM programs on demand. Packages like {admiral} become optional — a style preference rather than a requirement. Just describe what you need and let the model figure it out.

The benchmark data from pharma-skills tells a different story.

@online{dickinson2026,
  author = {Dickinson, Jeff},
  title = {Why We Still Need \{Admiral\} in an Age of {AI}},
  date = {2026-06-14},
  url = {https://pharmaverse.github.io/blog/posts/2026-06-14_admiral_in_age_of_ai/admiral_in_age_of_ai.html},
  langid = {en}
}